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An  analysis  of  the  stability  properties  of  uncertain  multi variable  control  systems 
in  the  frequency  domain  is  presented.  Necessary  and  sufficient  stability  criteria  are 
reviewed  along  with  singular  value  scaling  techniques  for  characterizing  permissible 
uncertainties.  Such  scaling  methods  have  become  widely  accepted  tools  for  the  anal- 
ysis of  control  systems  in  the  presence  of  structured  uncertainties.  Included  in  this 
study  are  the  general  block  similarity  scaling  techniques  advanced  by  Doyle  and  the 
nonsimilarity  scaling  approach  of  Kouvaritakis  and  Latchman. 

For  element-by-element  structured  uncertainties,  both  scaling  methods  reliably 
compute  Doyle's  structured  singular  value,  ft,  which  provides  an  indication  of  sys- 
tem stability.  However,  the  similarity  scaling  formulation  has  the  disadvantage  of 
expanding  annxn  system  matrix  to  an  n2  x  n2  matrix  requiring  n2  -  1  optimization 
variables  to  compute  fi.  Using  nonsimilarity  scaling,  the  system  size  remains  n  with 
the  additional  benefit  of  requiring  only  2(rc  -  1)  optimization  variables. 
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The  results  of  this  work  show  that  for  scalar  uncertainties,  the  structure  may  be 
exploited  to  yield  a  similarity  scaling  method  which  requires  no  more  than  the  2(n-l) 
optimization  variables  needed  for  nonsimilarity  scaling.  Substantial  savings  in  floating 
point  operations  are  observed  for  various  system  sizes  enhancing  the  capability  of  this 
method  for  analysis  and  iterative  design.  A  similar  reduction  in  optimization  variables 
is  shown  to  hold  for  the  important  class  of  general  block  structured  uncertainties.  This 
reduction  leads  to  a  complete  solution  for  the  2  x  2  block  uncertainty  problem. 

In  addition,  a  direct  relationship  between  the  similarity  and  nonsimilarity  scaling 
matrices  is  presented.  This  direct  relationship,  coupled  with  a  reduction  of  opti- 
mization variables  like  that  shown  for  similarity  scaling,  provides  a  more  efficient 
implementation  of  the  Fan  and  Tits  vector  optimization  method  for  computing  fi. 

For  systems  with  repeated  maximum  singular  values,  both  similarity  and  nonsim- 
ilarity scaling  procedures  generally  fail  to  calculate  /j,  exactly.  By  invoking  the  major 
principal  direction  alignment  principle  of  Kouvaritakis  and  Latchman,  a  2(q  -  1)- 
dimensional  optimization  problem  is  proposed  for  estimating  ft  where  q  represents 
the  multiplicity  of  the  maximum  singular  value.  In  all  numerical  experience  to  date, 
these  estimates  and  the  corresponding  exact  values  of  pt  have  agreed  within  three 
significant  figures.  Application  of  these  techniques  to  design  is  discussed. 
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CHAPTER  1 
INTRODUCTION  TO  SYSTEM  UNCERTAINTY 

1.1     Introduction 

The  first  task  in  any  control  design  is  the  modelling  process  through  which  a  math- 
ematical representation  of  the  system  is  developed.  This  model  may  be  constructed 
theoretically  through  some  knowledge  of  the  physical  laws  involved  or  empirically 
by  characterizing  experimental  data  collected  over  the  range  of  operating  conditions 
[1].   Theoretical  models  may  assume  linearity  around  some  nominal  operating  point 
while  neglecting  nonlinear  factors  that  may  dominate  away  from  this  nominal  point. 
Similarly,  the  range  of  operating  conditions  chosen  for  the  empirical  model  may  not 
be  sufficiently  exhaustive  to  adequately  characterize  system  behavior.  For  both  mod- 
elling approaches,  the  effects  of  component  aging,  temperature  and  pressure  varia- 
tions, manufacturing  tolerances  and  countless  other  unknown  factors  combine  to  cast 
doubt  on  the  accuracy  of  a  particular  system  model.  It  is  therefore  critical  to  ensure 
that  a  controller  designed  for  a  nominal  model  does  not  lose  required  stability  and 
performance  properties  when  applied  to  the  real- world  system    [2]. 

For  Single- Input  Single-Output  (SISO)  systems,  model  uncertainty  problems  have 
typically  been  addressed  by  ensuring  that  adequate  gain  and  phase  margins  exist 
throughout  the  range  of  operating  conditions.  These  margins  could  then  absorb 
the  detrimental  effects  of  the  uncertainties  without  sacrificing  stability  requirements. 
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For  Multi-Input  Multi-Output  (MIMO)  systems,  the  classical  definitions  of  gain  and 
phase  margins  do  not  apply  because  of  complex  system  interactions  [3].  Examples  of 
multivariate  systems  include  chemical  plants,  nuclear  facilities  and  high  performance 
aircraft  where  unstable  operation  may  result  in  the  loss  of  expensive  equipment  or 

even  lives    [4,  5]. 

y 

This  dissertation  reviews  some  of  the  current  frequency  domain  techniques  em- 
ployed in  the  analysis  of  multivariable  systems  with  uncertain  plant  models.  Of  par- 
ticular importance  are  those  approaches  that  address  structured  uncertainties  because 
they  generally  produce  a  much  less  conservative  stability  analysis  than  those  dealing 
strictly  with  unstructured  uncertainties.  Alternative  formulations  of  these  approaches 
are  developed  offering  substantial  improvements  in  computational  efficiency. 

Following  this  introduction  is  a  brief  discussion  of  the  historical  treatment  of  un- 
certainties for  feedback  systems  as  well  as  a  list  of  standard  notation  used  throughout 
this  work.  Chapter  2  summarizes  the  classical  Nyquist  stability  criterion  for  SISO 
systems  as  well  as  the  generalized  nyquist  criterion  for  MIMO  systems.  Singular  value 
techniques  that  provide  sufficient  stability  criteria  are  reviewed  in  Chapter  3  while 
Chapter  4  covers  various  scaling  techniques  that  recover  the  necessity  of  the  stability 
criteria  under  most  conditions.  The  main  contributions  of  this  dissertation  appear 
in  Chapters  5,  6,  7,  and  8  along  with  several  examples  that  illustrate  the  advantages 
offered  by  these  new  results.  Finally,  Chapter  9  summarizes  the  work  and  introduces 
some  of  the  promising  directions  for  future  robust  control  research. 
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1.2     Historical  Treatment  of  Uncertainty 

One  of  the  first  efforts  to  account  for  model  uncertainties  during  control  system 
design  was  by  H.  S.  Black  of  Bell  Laboratories  in  1927  [6].  Black  introduced  the 
concept  of  feedback  to  eliminate  amplifier  distortion  in  long  distance  telephone  com- 
munications. Amplifiers  of  that  time  required  hourly  adjustments  of  bias  currents 
resulting  in  unacceptable  manpower  demands  [7].  Black's  feedback  amplifier  proved 
to  be  almost  completely  immune  to  amplifier  uncertainties  caused  by  nonlinearities 
and  temperature  and  aging  changes.  One  drawback  to  the  new  amplifier  design  was 
an  occasional  self-oscillation  or  "singing"  at  certain  loop  gain  settings.  In  response 
to  this  sometimes  severe  problem,  another  Bell  Laboratories  scientist,  H.  Nyquist, 
developed  the  now  famous  Nyquist  stability  criterion  relating  closed  loop  system  sta- 
bility to  open  loop  frequency  response  information  [8].  From  this  stability  theory 
came  the  pioneering  work  of  H.  W.  Bode  concerning  the  issue  of  stability  robustness 
in  controller  design  [9].  The  resulting  concepts  of  gain  and  phase  margins  for  SISO 
systems  provided  the  design  engineer  with  an  measure  of  stability  in  the  presence  of 
uncertainties. 

The  Nyquist  stability  criterion  is  recognized  as  a  major  advance  in  the  design  of 
stable  SISO  control  systems  and  it  forms  the  basis  of  classical  control  theory.  How- 
ever, events  of  the  1950s  and  1960s  shifted  the  emphasis  of  control  theory  research 
from  the  frequency  domain  to  the  time  domain  representation  as  worldwide  atten- 
tion focused  on  manned  and  unmanned  rocket  guidance  and  control.  With  almost 
unlimited  budgets  and  well  defined  models,  the  problems  associated  with  system  un- 
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certainty  became  (at  least  temporarily)  less  important.  However,  the  late  1970s  saw 
renewed  interest  in  frequency  domain  analysis  and  design  as  the  precise  plant  models 
required  for  optimal  state-space  methods  simply  could  not  be  determined  for  many 
important  systems  due  to  numerous  plant  uncertainties    [10]. 

Finally,  in  1977  a  multivariable  analogue  to  the  Nyquist  stability  criterion  known 
as  the  characteristic  loci  method  (CLM)  was  developed  [11].  Introduced  by  Mac- 
Farlane  and  Postlethwaite,  the  CLM  utilizes  frequency  dependent  plots  of  transfer 
matrix  eigenvalues  (called  loci)  to  characterize  MIMO  system  stability.  In  its  original 
formulation,  the  CLM  does  not  provide  easily  discernable  information  on  stability 
margins  because  slight  perturbations  in  the  transfer  matrix  can  cause  large  shifts  in 
the  loci.  This  difficulty  led  to  the  use  of  singular  value  techniques  as  a  means  of 
extending  the  CLM  to  account  for  system  uncertainty. 

For  unstructured  uncertainties,  which  presume  no  knowledge  of  an  uncertainty 
other  than  a  norm  bound,  Doyle  and  Stein  have  shown  that  singular  value  bounds 
can  provide  necessary  and  sufficient  stability  conditions  for  uncertain  multivariable 
systems  [12].  Unfortunately,  this  formulation  fails  to  take  advantage  of  any  available 
knowledge  concerning  the  possibly  well-defined  structure  of  the  uncertainty.  The 
controller  must  then  accommodate  uncertainties  that  may  be  physically  impossible 
leading  to  an  overly  conservative  design    [13]. 

Safonov  applied  the  concept  of  similarity  scaling  to  singular  value  analysis  for 
the  case  of  diagonally  structured  scalar  uncertainties  [14].  This  scaling  idea  was 
then  extended  by  Doyle  to  consider  the  important  class  of  general  block  structured 


perturbations  which  allow  norm  bounded  uncertainty  blocks  to  be  of  arbitrary  dimen- 
sion [15].  Doyle's  concept  of  the  "structured  singular  value"  (denoted  ft)  as  a  means 
of  determining  the  set  of  permissible  structured  uncertainties  for  system  stability  has 
become  a  widely  accepted  tool  in  the  design  of  robust,  multivariate  control  systems. 
Additional  scaling  techniques  such  as  the  nonsimilarity  approach  of  Kouvaritakis  and 
Latchman  have  evolved  to  address  element-  by-element  structured  uncertainties  and, 
in  fact,  both  similarity  and  nonsimilarity  scaling  methods  have  been  shown  to  pro- 
duce necessary  and  sufficient  stability  conditions  in  terms  of  fi  for  most  systems  with 
structured  uncertainties    [16,  17]. 

These  singular  value  techniques  have  recently  been  extended  to  the  areas  of  #°° 
and  //-Synthesis:  both  of  which  allow  a  great  deal  of  flexibility  in  the  satisfaction  of 
controller  design  requirements  while  retaining  stability  and  performance  robustness 
properties  [18,  19].  The  main  emphasis  of  this  work  will  therefore  be  directed  to- 
wards the  singular  value  analysis  techniques  including  new  results  that  enhance  their 
application  in  the  analysis  of  robust  multivariable  control  systems. 

1.3     Notation 

The  following  notational  convention  will  apply  unless  otherwise  stated. 
CnXm    :     The  set  of  complex  matrices  with  n  rows  and  m  columns. 
sjjnxm    .     q^g  get  Q£  reaj  matrjceg  w^n  n  rows  anci  m  columns. 

j  :     The  square  root  of  —1. 

Im(a)  :     The  imaginary  portion  of  complex  element  a. 

Re(a)   :     The  real  portion  of  complex  element  a. 


arg(a)  :  The  argument  of  complex  element  a. 

\a\         :  The  absolute  value  of  element  a. 

a  :  The  complex  conjugate  of  scalar  a. 

AH        :  The  complex  conjugate  transpose  of  matrix  A. 

\\A\\P    :  The  p-novm  of  matrix  A  (p  =  2  unless  noted  otherwise). 

\\A\\f    ■  The  Frobenius  norm  of  matrix  A. 

A+        :  Matrix  A  with  elements  replaced  with  their  absolute  values. 

A-1       :  Inverse  of  matrix  A. 

det{A}:  The  determinant  of  matrix  A. 

\i(A)    :  The  zth  eigenvalue  of  matrix  A. 

p(A)     :  The  spectral  radius  of  A, 

(Ji{A)    :  The  iih.  singular  value  in  magnitude. 

a(A)     :  The  singular  value  with  maximum  magnitude  (a  =  ax). 

fi(A)     :  The  structured  singular  value  of  matrix  A. 

V         :  The  family  of  diagonal  matrices  with  positive,  real  entries. 

U  :  The  family  of  diagonal  unitary  matrices. 

inf        :  Infimum. 

sup       :  Supremum. 

max      :  Maximum. 

■  :  The  completion  of  a  proof  or  discussion. 

CPU     :  Central  processing  unit. 


CHAPTER  2 
FREQUENCY  DOMAIN  ANALYSIS 

2.1     Transfer  Matrix  Representation 

Any  real  world  system  can  be  characterized  by  a  relationship  between  system 
inputs  and  outputs.  Given  a  system  with  output  y(s)  and  input  u(s),  the  transfer 
function  G(s)  relates  the  two  in  the  manner 

y(s)  =  G(s)u(s), 

For  nonscalar  G(s),  the  off-diagonal  elements  of  G(s)  produce  the  system  interac- 
tion that  complicates  multivariable  control  systems.  The  identification  process  to 
determine  a  transfer  matrix  begins  by  injecting  known  inputs  u(s)  into  the  plant  and 
measuring  the  resulting  output  y(s).  An  alternative  approach  begins  with  the  time 
domain  state  space  equations  usually  derived  through  some  knowledge  of  the  plant's 
physics.  The  state  space  representation  may  be  expressed  as 

x(t)  =  Ax(t)  +  Bu{t),       y(t)  =  Cx(t)  +  Du(t) 


where  A,  B,  C,  and  D  are  real,  possibly  time  varying  matrices.  Through  the  Lapl 
transform,  the  unique  transfer  matrix  representation  may  then  be  written 


ace 


as 


G(s)  =  C(sI-A)-1B  +  D. 


Figure  2.1:     State  Space  System  Representation. 


Figure  2.2:     Transfer  Matrix  System  Representation. 


2-2     Characteristic  Loci  and  the  Generalized  Nyquist  Criterion 

Two  assumptions  are  made  at  this  point  concerning  the  transfer  matrices  used 
throughout  this  work.  The  first  assumption  is  that  G(s)  is  square  and  rational. 
The  second  assumption  is  that  G(s)  contains  no  hidden  unstable  modes.  Also,  only 
unity  feedback  will  be  addressed  explicitly  although  any  uncompensated  G(s)  may 
simply  be  rewritten  as  G(s)  =  K(s)G(s)  where  K{s)  is  some  precompensator.  The 
formulation  would  then  continue  with  G(s)  instead  of  G(s). 

Denote  g(s)  =  ^  as  the  open-loop  transfer  function  of  a  scalar  system  where 
n(s)  and  d(s)  represent  the  numerator  and  denominator  polynomials  respectively. 
Then  g(s)  may  be  related  to  the  return  difference  transfer  function,  /($),  by 

ft    \         i     ,       t    ^         n(s)  +  d(s) 

f(s)  =  l+g(s)=         >—±  >.  (2.1) 

The  term  n{s)  +  d(s)  from  Equation  2.1  actually  corresponds  to  the  closed-loop 
pole  polynomial,  pc(s),  under  unity  feedback  while  d(s)  corresponds  to  the  open-loop 
pole  polynomial,  p0,  of  g(S).  Therefore  f(s)  may  be  rewritten  as  the  ratio  of  the 
closed-loop  to  open-loop  pole  polynomials  or 

m  =  l+9(s)  =  ^f\.  (2.2) 

Application  of  the  principle  of  the  argument  allows  the  development  of  a  graphical 
stability  test  for  Equation  2.2.  Denote  the  open  and  closed-loop  right-half  plane  poles 
as  n0  and  nc  respectively.  Mapping  g{s)  onto  the  classical  Nyquist  D-contour  requires 
that  the  number  of  clockwise  critical  point  (-1  +  jO)  encirclements  equals  nc  -  n0. 
Closed-loop  stability  implies  that  nc  =  0  so  the  number  of  critical  point  encirclements 
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for  a  stable  system  must  be  -n0.  This  result  forms  the  basis  for  the  classical  Nyquist 
stability  criterion. 

Extending  these  results  to  the  multivariable  case,  MacFarlane  and  Postlethwaite 
generalized  Equation  2.2  to  the  form 

det{F(s)}  =  det{[I+  G(s)}}  =  a^\ 

Po{s) 

where  F(s)  is  the  multivariable  return  difference  matrix  and  a  is  a  scalar  constant 
[11].  Rather  than  mapping  g(s)  directly  as  in  the  SISO  case,  the  characteristic  loci, 
defined  as  the  Nyquist  D-contour  map  of  the  frequency  dependent  eigenfunctions, 
gi{s),  are  mapped.  These  eigenfunctions  are  solutions  to  the  equation 

det{[9i(s)I  -  G(s)]}  =  0 

and  hence  form  algebraic  functions  of  the  elements  of  G(s).  Also,  the  9i(s)  are  analytic 
everywhere  except  where  two  or  more  of  the  g{(s)  coincide.  At  such  points,  the 
separate  layers  of  the  Riemann  surface  which  forms  the  domain  of  the  gi(s)  are  simply 
joined  together  forming  closed  curves.  Once  again  the  principle  of  the  argument  may 
be  applied  to  these  closed  curves  to  form  the  basis  of  the  generalized  Nyquist  criterion 
stated  in  Theorem  1.1    [11]. 

Theorem  1.1:  A  multivariable  system  with  open-loop  transfer  matrix  G(s)  con- 
taining only  controllable  and  observable  modes  will  be  stable  under  unity  feedback  if 
and  only  if  the  net  sum  of  counter-clockwise  characteristic  loci  encirclements  of  the 
critical  point  (-1  +  j0)  is  equal  to  the  number  of  right-half  plane  open-loop  poles  of 
G(s). 
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This  important  theorem  allows  for  a  complete  graphical  stability  test  of  multi- 
variable  systems  in  the  frequency  domain  where  no  model  uncertainty  is  present. 
Unfortunately,  SISO  concepts  of  system  robustness  based  on  gain  and  phase  mar- 
gins do  not  extend  to  the  multivariable  case  as  the  eigenfunctions  may  display  large 
shifts  from  relatively  small  perturbations  in  the  elements  of  G(s).  The  next  section 
addresses  the  stability  of  systems  in  the  presence  of  uncertainties. 

2.3      Necessary  and  Sufficient  Stability  Conditions 

The  generalized  Nyquist  criterion,  developed  from  the  characteristic  loci  method, 
provides  both  necessary  and  sufficient  stability  conditions  for  the  nominal  plant  de- 
noted G0(s)  where  a  =  jlo  on  the  Nyquist  D  contour.  Stability  is  guaranteed  if  and 
only  if  the  number  of  counter-clockwise  critical  point  encirclements  by  the  character- 
istic loci  equals  the  number  of  unstable  open  loop  poles.  Assuming  that  the  nominal 
plant  itself  is  stable  (either  alone  or  following  some  type  of  control  implementation), 
instability  in  the  perturbed  plant,  Gp(s),  implies  a  change  in  the  number  of  critical 
point  encirclements.  If  G0(s)  and  Gp(s)  have  the  same  number  of  of  open-loop  unsta- 
ble poles,  then  for  such  a  change  to  occur,  at  least  one  eigenvalue  of  Gp(s)  must  pass 
through  the  critical  point  (-1  +  ;0).  Restricting  the  nominal  and  perturbed  plants 
to  the  same  number  of  open-loop  poles  requires  uncertainty  matrix  A  to  be  stable 
transfer  matrix.  While  this  requirement  does  limit  allowable  uncertainty  structures  to 
those  with  no  poles  in  the  right-half  plane,  Foo  has  shown  that  unstable  perturbations 
can  generally  be  decomposed  into  two  stable  perturbations  allowing  the  analysis  to 
continue  with  a  possible  increase  in  conservatism  [20]. 
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For  additive  uncertainties,  the  necessary  and  sufficient  stability  criterion  for  an 
uncertain  plant  Gp(s)  =  (G0(s)+A(s)),  where  A(s)  represents  a  frequency  dependent 
uncertainty  matrix,  becomes 

KGP(s))  #  -1. 
Substituting  back  the  nominal  plant  and  A(s),  the  criterion  may  be  written  as 

\(G0{s)  +  A(s))  ^ -1 


or 


X(G0(s)  +  A(s)  +  I)  ^  0.  (2-3) 

For  any  matrix  A,  the  determinant  of  A  equals  the  product  of  the  eigenvalues  of  A 
so  Equation  2.3  becomes 

det{(G0(s)  +  A(s)  +  I)}  ^  0 
or 

det{(G0(s)  +  1}  ■  det{(I  +  Gois))-1  A(s)  +  1}  ±  0.  (2.4) 

The  nominal  plant  is  known  to  be  stable  and  hence  cannot  have  any  eigenvalues 
passing  through  the  critical  point.  Applying  this  fact,  Equation  2.4  may  be  simplified 


to  the  form 


X(I  +  M(s)A(s))  ^  0        where  M(s)  =  (I  +  G^s))'1 


or 


\(M(s)A(s))  jL  -1.  (2-5) 
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Typically,  some  form  of  norm  bound  on  either  the  entire  matrix  (unstructured  uncer- 
tainties) or  individual  elements  (structured  uncertainties)  can  be  estimated,  leaving 
the  phase  information  to  vary  freely.  This  freedom  in  choosing  the  phase  of  A  allows 
Equation  2.5  to  be  written  as  the  necessary  and  sufficient  stability  condition 

Supp(M(s)A(s))<l  (2.6) 

where  p  denotes  the  spectral  radius.  It  should  be  noted  that  while  the  preceding 
argument  addresses  additive  uncertainties,  a  similar  approach  can  be  applied  to  mul- 
tiplicative uncertainties  where  Gp(s)  =  G0(s)(I  +  A(s)).  Multiplying  on  the  right  by 
G0(s)  gives 

Gp(s)  =  GD(s)  +  G0(s)A(s) 
or 

Gp(s)  =  G0(s)  +  A(s)     where  A(s)  =  G0(s)A(s). 

Although  Equation  2.6  provides  necessary  and  sufficient  stability  conditions  for 
uncertain  systems,  it  is  quite  difficult  to  compute  because  the  entire  range  of  A  must 
be  considered.  Furthermore,  the  solution  lacks  convexity  allowing  for  the  possibility 
of  multiple  maxima.  This  unfortunate  situation  greatly  complicates  attempts  to  find 
the  true  supremum  of  p(M(s)A(s))  over  A(s).  Because  of  the  difficulties  involved 
in  solving  Equation  2.6,  singular  value  techniques  were  developed  to  provide  upper 
bounds  on  permissible  model  uncertainties.  Before  introducing  these  techniques,  a 
mathematical  description  of  several  important  uncertainty  classes  will  be  presented. 
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2.4     Uncertainty  Classes 

So  far,  the  frequency  dependent  uncertainty  matrix,  A(s),  has  been  used  to 
represent  the  additive  or  multiplicative  uncertainty  present  in  control  system.  No  in- 
formation has  been  provided  on  the  structure  of  A(s)  which  turns  out  to  be  quite  im- 
portant in  the  actual  mathematical  analysis  of  system  stability.  Rather  than  propos- 
ing one  all-encompassing  model  for  A (s)  to  account  for  system  uncertainty,  three 
representations  will  be  presented  with  varying  degrees  of  information  required  for 
each  one.  In  developing  these  three  models,  two  main  objectives  are  observed: 

1.  The  model  should  handle  all  information  available  concerning  the  system's  ac- 
tual uncertainties.  All  impossible  uncertainty  structures  should  be  excluded  to 
reduce  conservatism  of  the  stability  analysis. 

2.  The  model  should  be  as  simple  as  possible  so  that  the  analysis  process  is  not 
unnecessarily  complex.  Simplified  models  also  encourage  an  interactive  design 
process  so  various  control  designs  may  be  quickly  generated  and  compared. 

Naturally,  these  objectives  produce  conflicting  requirements  and  careful  consideration 
must  be  used  to  select  an  appropriate  model.  In  any  event,  a  conservative  stability 
analysis  is  preferable  over  one  that  is  simple  but  erroneous. 

The  three  uncertainty  models  considered  here  include  unstructured  uncertainties, 
element-by  element  structured  uncertainties,  and  block  diagonal  structured  uncer- 
tainties. In  each  of  these  classes,  the  uncertainty  matrix  A(s)  is  complex  with  some 
form  of  norm  bound  placed  on  the  magnitude  of  the  matrix  or  matrix  elements.  The 
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simplest  uncertainty  class  to  describe  and  manipulate  concerns  unstructured  uncer- 
tainties which  may  be  mathematically  represented  as 

DU  =  {A\  ||A|ja  <  8  €  3?+} 
where  ||  ■  ||2  represents  the  standard  matrix  2-norm  and  8  a  real  scalar.  (Here,  and 
in  the  remainder  of  this  work,  frequency  dependence  will  be  implied).  From  the 
definition,  it  is  apparent  that  no  information  concerning  the  inner  structure  of  A  is 
accessible.  This  effectively  places  a  SISO  bound  on  a  MIMO  system  greatly  hinder- 
ing efforts  to  reduce  conservatism  if  information  about  the  inner  structure  of  A  is 
available. 

The  second  uncertainty  class,  element-by-element  structured  uncertainties,  pro- 
vides a  rich  source  of  information  about  A.  This  class  is  defined  as 

Ds  =  {A  =  {8lJ}\\8lj\<plJ  G&+} 

where  a  magnitude  bound  is  placed  on  each  element  of  uncertainty  matrix  A  allowing 
information  on  the  inner  structure  of  A  to  remain  a  part  of  the  stability  analysis 
process.  Class  Ds  provides  a  more  realistic  representation  of  real  world  uncertainty 
than  that  of  the  unstructured  uncertainty  class. 

To  illustrate  the  advantages  of  structured  uncertainties  over  unstructured  uncer- 
tainties, consider  some  process  plant  with  two  independent  control  valves  whose  flow 
rate  is  known  within  ±10%  [21].  As  a  structured  uncertainty,  A  could  be  correctly 
represented  as 

8X     0 


A  = 


0     8o 


\8i\  <  0.1. 
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Placing  a  norm  bound  of  ||  A||2  <  .1  as  required  for  unstructured  uncertainties  allows 
a  number  of  impossible  structures  like 


A  = 


0     0 
0.1    0 


r                        -I 

1 

,     or    A  =  -  • 

2 

0.1    0.1 
0.1    0.1 

Since  there  are  two  uncertainties,  not  one  or  four,  these  do  not  correspond  to  any  real- 
izable uncertainty  configuration  so  a  stability  analysis  accounting  for  such  impossible 
structures  would  be  unnecessarily  conservative. 

The  third  uncertainty  class,  block  structured  uncertainties,  is  actually  a  superset 
of  the  first  two  classes.  Mathematically  represented 


as 


Db  =  {A  |  HAtfHa  <  %  €  &+} 

the  class  Dh  may  have  unstructured  submatrices,  Ay,  with  arbitrary  dimension.  For 
the  case  where  size(Al3)  =  1x1,  this  class  reduces  simply  to  Ds.  Likewise,  for 
*tze(Ay)  =  n  x  n  where  n  is  the  size  of  A,  class  Db  reduces  to  class  Du.  Block 
structured  uncertainties  advanced  from  work  done  by  Doyle  [15]  and  they  have  become 
an  important  area  of  study  due  to  their  quite  general  nature. 

Most  real- world  system  uncertainties  may  be  cast  in  one  of  the  three  classes  dis- 
cussed above.  However,  choosing  the  particular  class  for  a  specific  problem  is  seldom 
straightforward  and  assigning  actual  magnitudes  for  the  elements  of  the  chosen  struc- 
ture can  be  more  difficult  still.  Comparing  the  behavior  of  nominal  plant  G0  with  an 
uncertainty  matrix  A  to  that  of  the  actual  system  provides  insight  into  the  accuracy 
of  the  chosen  model  although  a  large  number  of  candidate  plants  and  uncertainty 
structures  may  satisfactorily  match  experimental  data  gathered  during  the  identifi- 
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cation  process.   Narrowing  down  which  of  these  candidate  model  combinations  best 

approximates  the  actual  system  becomes  a  process  of  model  invalidation  as  those 
models  that  fail  to  satisfy  some  matching  criteria  are  eliminated.  No  discussion  of 
how  a  designer  should  choose  one  model  over  another  will  be  presented  here  since 
this  is  a  complete  field  unto  itself.  However,  a  description  of  the  problem  by  Smith 
and  Doyle  may  be  found  in    [22]. 

While  most  real- world  systems  may  be  characterized  by  one  or  more  of  the  uncer- 
tainty classes  described  here,  it  must  be  noted  that  other  models  offering  improved 
model  agreement  exist  for  some  specific  problems.  Most  notable  of  these  other  mod- 
els are  those  that  deal  with  combinations  of  both  real  and  complex  uncertainties  and 
those  that  address  repeated  or  dependent  uncertainty  structures.  In  the  three  uncer- 
tainty structures  discussed  above,  unstructured,  element-by-element  structured,  and 
general  block  structured,  some  form  of  bound  was  placed  on  the  magnitude  of  the 
uncertainty.    Such  magnitude  bounds  allow  the  phases  of  the  uncertainty  elements 
or  blocks  to  vary  freely  between  0  and  2tt  so  that  worst  case  situations  can  be  ad- 
dressed. For  many  uncertainty  sources  like  unmodelled  dynamics,  this  phase  freedom 
is  mandatory  to  establish  necessary  and  sufficient  stability  conditions.    However,  if 
some  or  all  elements  of  an  uncertainty  structure  are  known  to  be  strictly  real,  this  al- 
lowable phase  variation  may  provide  for  an  overly  conservative  stability  analysis  [23]. 
A  method  for  treating  strictly  real  uncertainties  was  proposed  by  De  Gaston  and 
Safonov  [24].  This  method  transforms  the  perturbed  plant  Gp  into  a  series  of  convex 
hulls  in  the  complex  plane  and,  through  an  iterative  process,  determines  a  stabil- 
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ity  margin  that  indicates  the  largest  scalar  multiplier  of  A  for  which  stability  is 

guaranteed.  This  approach  provides  necessary  and  sufficient  stability  conditions  but 

only  if  an  infinite  number  of  iterations  are  performed.    Fortunately,  usable  results 

are  generally  obtained  for  a  finite  number  of  iterations  though  the  method  is  still 

computationally  intensive. 

For  uncertainties  consisting  of  combinations  of  real  and  complex  elements,  alter- 
native formulations  of  the  structured  singular  value  have  been  proposed  by  several 
authors  [23,  25,  26].  While  these  approaches  appear  to  reduce  the  conservatism  over 
methods  that  consider  only  complex  uncertainty  structures,  a  number  of  theoretical 
and  computational  challenges  must  be  resolved  before  they  can  be  reliably  applied  to 
engineering  problems. 

In  addition  to  special  uncertainty  models  addressing  real  elements,  the  case  of 
repeated  or  dependent  complex  uncertainty  models  requires  special  consideration. 
This  class  is  actually  a  subset  of  the  other  structured  uncertainty  classes  with  the 
distinction  that  two  or  more  of  the  elements  or  blocks  of  A  are  constrained  such  that 
their  phases  must  always  be  the  same  [27].  Thus,  while  the  phases  may  vary  between 
0  and  27r,  they  may  not  do  so  independently  of  one  another. 

While  the  real  and  repeated  uncertainty  element  models  provide  an  improved 
stability  analysis  for  their  particular  cases,  the  two  structured  uncertainty  classes 
presented  earlier  account  for  a  wide  range  of  important  problems.  Therefore,  only 
their  treatment  will  be  considered  in  the  remainder  of  this  work. 


CHAPTER  3 
SINGULAR  VALUE  TECHNIQUES 

3-1       The  Singular  Value  Decomposition 

As  discussed  in  Section  2.3,  the  necessary  and  sufficient  stability  criterion 

sup^(ATA)<l  (3.1) 

is  nonconvex  making  it  difficult  if  not  impossible  to  determine  the  range  of  permissi- 
ble As.  Rather  than  giving  up  on  this  approach  for  uncertain  multivariable  systems, 
it  is  possible  to  continue  the  analysis  using  relationships  involving  the  singular  value 
decomposition.  These  relationships  start  out  as  conservative  upper  bounds  but  the 
results  of  Chapter  4  show  that  this  conservatism  can  generally  be  eliminated  com- 
pletely. The  singular  value  decomposition  is  defined  for  any  matrix  A  £  CnXn 


as 


A  =  XVYH  (3.2) 

where  X  and  Y  are  unitary  matrices  the  columns  of  which  form  the  left  and  right 
singular  vectors  of  A  respectively.  The  diagonal  matrix  £  contains  the  singular  values 
of  A  in  decreasing  order  of  magnitude.  The  singular  values  and  vectors  of  A  can  be 
determined  in  terms  of  the  eigenvalues  and  eigenvectors  of  the  hermetian  forms  AH A 
and  AAH  as 

AHAYi    =    <rfYi 

AAHXt    =    afXi.  (3.3) 
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Other  useful  relationships  of  the  singular  values  and  corresponding  singular  vectors 
include: 

AYi  =  diXi 

Y»AH  =  criAf 

AHXi  =  aiY 

XfA  =  a^/1.  (3.4) 

Definition  3.1  The  2-norm  of  A  €  CnXr\  ||A||2,  is  defined  as 

ll^2  =  oS»ir^  (3-5) 

o^xeCn    \\x\\2 
where  ||x||2  is  the  vector  2-norm  defined  as 

Also,  from  Definition  3.1  comes  the  important  inequality  equation 

\\Az\U  £  \\Ah\M2.  (3.6) 

The  following  theorem  relates  the  2-norm  of  a  matrix  to  its  maximum  singular  value. 
Theorem  3.1:  For  any  matrix  A  G  CnXn 

nA)  =  \\A\\2. 

Proof:  From  Equation  3.2,  matrix  A  may  be  written  as  A  =  XY,YH '.  The  2-norm 
of  a  matrix  is  invariant  to  unitary  transformations  so 

||A||2  =  \\XTYH\\%  =  HUlk 
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From  Definition  3.1,  this  may  be  written  as 

ii  4  ii       urn  \\^xh  \lori\xi\2  +  ---  +  °l\xn\2 

\\A-h  =   p  2  =    max     '     „"     =    max   - — . 

o^ec-   ||aj||2       o^eo      ^2  +  . . .  =  j^jj 

Choosing  x  =  ei,  the  first  standard  basis  vector  gives  lBjjjjjb  =  ffi(A)  while  any  other 
choice  for  x  gives  J»  <  ^(A).  Therefore,  ||A||2  =  a^A)  =  7(A).  ■ 

Using  this  result,  the  following  lemma  relates  a  (A)  to  p(A). 

Lemma  3.1  For  any  matrices  A,  B  G  CnXn, 

/>(A5)<<f(A£).  (3.7) 

Proof:  From  the  eigenvalue/eigenvector  equation  we  have 

ABWi  =  XiWi. 

where  Aj  is  an  eigenvalue  of  (AB)  and  W{  the  corresponding  eigenvector.  Taking  the 
2-norm  of  both  sides  gives 

||A5^||2  =  ||A!^||2  =  |Ai|||^||2. 
Equation  3.6  allows  this  to  be  written  as  an  inequality 

\\AB\\2\\Wl\\2>\\ABWl\\2  =  \Xi\\\Wl\\2. 
Applying  Theorem  3.1  along  with  the  definition  of  the  spectral  radius  gives 

a(AB)  >  p(AB).    m 


Finally,  the  2-norm  of  a  matrix  product  is  related  to  the  product  of  the  2- 
as  shown  in  the  following  lemma. 


norms 


Lemma  3.2  For  A,B  £  CnXn 

a(AB)<a(A)a(B).  (3J 

Proof:  Again,  from  Definition  3.1 

UBh=    max   M*<    max   MMMl  ,  M>   max   M^ 
<¥*ec«     ||x||2         oA6C»        ||x||2  l|/1||2o™tc»    IUIL 
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or 


P5||3  <  ||A||2||5||a.    - 

Replacing  A  and  5  with  M  and  A  respectively  in  Equations  3.7  and  3.8,  the  following 
expression  must  hold: 

p(MA)  <  a(MA)  <  a(M)a(A). 

From  these  relationships,  the  necessary  and  sufficient  stability  criterion  of  Equation 
3.1  can  be  written  as  the  following  sufficient  stability  conditions 

supa(M)a(A)  <  1  (3_9) 

and 

supo^AfA)  <  1.  (3  10) 

Since  Equations  3.9  and  3.10  only  guarantee  sufficient  stability  conditions  for 
structured  uncertainties,  additional  manipulations  are  required  to  reduce  conser- 
vatism with  the  intent  of  regaining  the  necessary  stability  condition.  Two  techniques 
that  address  this  problem  are  nonsimilarity  scaling  introduced  by  Kouvaritakis  and 
Latchman  [16]  and  similarity  scaling  advanced  by  Doyle  [15].  Both  techniques  rely 
on  the  fact  that  while  the  eigenvalues  (and  hence  the  spectral  radius)  of  a  matrix 
unaffected  by  similarity  transformations,  the  singular  values  are  not  necessarily 


are 
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preserved  and  may,  in  fact,  be  reduced.  The  nonsimilarity  scaling  technique  will  be 
discussed  first  with  the  analysis  limited  to  element  by  element  bounded  uncertainties. 

3-2     Structured  Uncertainties  with  Nonsimilarity  Scaling 

The  nonsimilarity  scaling  formulation  starts  with  the  necessary  and  sufficient 
stability  criterion  of  Equation  3.1.  Introducing  positive,  diagonal  matrices  R  and  L 
to  this  expression  gives  the  equivalent  stability  condition 

supp(MA)  =  sup  pfR-1  ML'1  LAR)  <  1. 

A  A 

This  may  be  rewritten  in  the  form  of  Equation  3.9  giving  the  sufficient  stability 
criterion 

supcr(ir1MJL-1)<f(Z,A.R)  <  1.  (3.11) 

A 

Equation  3.11  no  longer  conforms  to  the  similarity  transformation  structure:  hence 
the  name  "nonsimilarity  scaling."  However,  the  free  elements  of  R  and  L  still  allow 
the  conservatism  gap  between  the  maximum  singular  value  upper  bound  and  the 
spectral  radius  lower  bound  to  be  reduced.  Also,  the  choice  of  positive  values  for 
the  elements  of  R  and  L  allows  a  simplification  based  on  the  following  lemmas  and 
Theorem  3.2. 

Lemma  3.3:  For  all  A  €  CnXn  and  B  £  $n*n  where  b{j  >  0  and  \al3\  <  %  for  all 
i,j,  with  x+  defined  as  (JxjJ,  •  •  • ,  |:r„|)T, 

||Ax||a  <  \\Bx+\\2  for  all  x  e  Cn.  (3.12) 
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Proof:  From  the  definition  of  a  vector  2-norm 


||Ar||2    =     \/\a11x1  +  ---  +  alnxn\2  +  ---  +  \anlx1  +  ---  +  annxn\2 
WBx+h    =     V7!  [&nW  +  ---  +  &in|sn|]  I2  +  ••■  +  !  [bni^l  +  ■  ■  ■  +  bnn\xn\]  |2. 

The  second  equation  contains  only  nonnegative  numbers  so  the  sum  of  the  individual 
terms  cannot  decrease.   The  first  equation,  however,  contains  elements  that  may  be 
negative  so  the  sum  of  the  individual  terms  may  decrease.   Since  the  elements  of  B 
are  greater  than  or  equal  to  the  corresponding  elements  of  A,  the  lemma  must 
hold.  ■ 

Lemma  3.4  Using  the  definition  of  x+  from  Lemma  3.3, 

\\x\\2  =  ||x+||2   for  all  x  €  Cn.  (3.13) 

Proof:  Again,  from  the  definition  of  the  vector  2-norm 

IMIa  =  \J\xi\2  +  ---  +  W2  =  \\x+b-    ■ 
Lemma  3.5  For  A  and  B  defined  as  in  Lemma  3.3, 

^-^  <  ||5||2  for  all  0  ±  x  €  Cn.  (3.14) 

II*||3  y  J 

■ 

Proof:  From  Definition  3.1, 

II  nil  \\Bx\\2 

\\B\\-2  =    max    "„    „"   . 

0#*6C»     ||x||2 

Since  all  elements  of  B  are  nonnegative,  this  may  be  combined  with  Lemmas  3.3  and 
3.4  to  give 

11511,=  max  J!^>M*.  m 

«¥*+6C7»    ||x+||2    -    \\x  U 
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Theorem  3.2:  Given  A   6  C"*»,  positive,  diagonal  matrices  L,R   e  R—  and 

P  e  »nx"  where  ^  >  0  and  |%|  <  Pij  for  all  ij, 

a(LAR)  <a(LPR).  (3J5) 

Proof:  This  proof  follows  directly  from  Lemmas  3.3,  3.4,  and  3.5  with  A  replaced 
by  LAR  and  B  replaced  by  LPR  so  that 

II  r  A  on  l|Aa;|U  \\LPRx+\\2 

\\Li\K\\2  =    max    '„    ,,"    <     max ±-li£  —     r,pmi       ■ 

<¥*eO    ||x||2    -o^+ec»       ||x+||2       ~  H^^lla-    ■ 

Theorem  3.2  and  the  definition  of  element-by-element  structured  uncertainties 
from  Section  2.4  allow  the  direct  substitution  of  P  for  A  in  Equation  3.11. 

Substituting  P  for  A  would  appear  to  add  additional  conservatism  to  Equation 
3.11.  However,  Lemma  3.6  shows  that  this  is  not  the  case. 

Lemma  3-6  Given  LAR  and  LPR  defined  as  in  Theorem  3.2,  the  following  equality 
must  hold: 

supw{LAR)  =  a(LPR).  (3>16) 

Proof:  The  freedom  of  the  phases  of  each  element  of  A  allows  StJ  =  |<%|  =  Pij 
for  all  i  J.  Therefore,  by  individually  adjusting  the  phases  of  each  element  of  A,  the 
equality  of  Equation  3.16  is  established.  ■ 

Since  P  always  represents  a  realizable  variation  of  A,  there  is  actually  no  added 
conservatism  by  substituting  P  for  A.  Therefore  the  sufficient  stability  condition  of 
Equation  3.11  becomes 

iniRa{R-lML-l)a{LPR)  <  1  (3.17) 
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over  all  design  frequencies.    This  infimization  over  L  and  R  requires  only  2(n  -  1) 

optimization  variables  because  one  diagonal  element  of  both  L  and  R  may  be  held 
constant  without  loss  of  generality.  For  the  case  of  distinct  a(R~lML-1),  Equation 
3.17  has  been  shown  to  converge  to  the  necessary  and  sufficient  condition  of  Equation 
3.1.  This  optimal  condition  is  characterized  by  the  major  principal  direction  align- 
ment (MPDA)  theorem  proposed  by  Kouvaritakis  and  Latchman  [17].  Due  to  the 
importance  of  this  theorem  it  will  be  discussed  in  detail  in  Chapter  4. 

3-3     Structured  Uncertainties  with  Similarity  Scaling 


les 


The  previous  section  outlined  the  nonsimilarity  scaling  technique  which  appl 
to  systems  with  element-by-element  structured  uncertainties.  A  second  approach  to 
solving  Equation  3.1  through  singular  value  techniques  involves  diagonal  similarity 
scaling  where  only  one  scaling  matrix  is  used.  Advanced  by  Doyle  [15],  this  technique 
requires  that  the  uncertainty  matrix  A  be  diagonal  with  norm  bounds  on  the  diagonal 
elements  or  blocks.  An  important  advantage  of  diagonal  similarity  scaling  is  its 
ability  to  address  the  general  block  structured  uncertainties  discussed  in  Section  2.4. 
While  the  requirement  for  diagonalized  uncertainties  appears  to  severely  restrict  the 
applicability  of  this  technique,  it  is  possible  to  diagonalize  any  uncertainty  matrix 
through  the  use  of  simple,  eigenvalue  preserving  transformations.  Properties  of  these 
transformations  will  be  discussed  in  Section  3.5. 

For  diagonal  uncertainty  matrix  Ad,  diagonal  scaling  matrix  D  may  be  introduced 
allowing  the  necessary  and  sufficient  stability  conditions  of  Equation  3.1  to  be  written 
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as 

supp(MAd)  =  sup^DMA^-1)  <  1  (3.I8) 

and  the  sufficient  condition  of  Equation  3.10  as 

supapMA^Zr1)  <  1.  (3.19) 

Taking  advantage  of  the  diagonal  nature  of  Ad  permits  a  decomposition  of  the  form 

A,  =  PdUd 

where  pdu  >  \8dii\  for  all  i  and  Ud  is  a  diagonal  unitary  matrix.  This  allows  Pd  to 
retain  the  magnitude  information  and  Ud  the  phase  information  of  Ad.  Substituting 
PdUd  for  Ad  in  Equation  3.19  gives  the  sufficient  stability  condition 

supp(MAd)  <  sup  cr( DM PdUdD'1)  <  1. 

Ad  ud 

Noting  that  diagonal  matrices  commute  with  each  other,  the  positions  of  Ud  and  D'1 
may  be  exchanged  making  the  condition 

supp(MAd)  <  sup  uiDMPdD^Ud)  <  1.  (3.20) 

Ad  Ud  y  J 

However,  the  2- norm,  and  hence  the  maximum  singular  value,  of  a  matrix  is  invariant 
to  unitary  transformations  removing  the  dependence  on  Ud.  Combining  M  and  Pd 
into  Ma  =  MPd,  the  condition  of  Equation  3.20  becomes 

sup  p{MAd)  <  MaiDMaD'1 )  <  1 .  (3.21 ) 

Ad  D 

Since  the  spectral  radius  of  MaUd  is  always  bounded  from  above  by  the  maximum 
singular  value  of  DMaD~l,  the  free  elements  of  D  may  be  used  to  reduce  the  conser- 
vatism gap  between  the  maximum  singular  value  upper  bound  and  the  spectral  radius 
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lower  bound  of  Equation  3.21.  For  the  general  case  where  a  stationary  point  occurs 
at  the  "inf,"  the  resulting  optimal  value  of  a{DMaD~l)  is  known  as  the  "structured 
singular  value"  of  Ma  (denoted  p(Ma)). 

The  original  definition  of  the  structured  singular  value  had  the  form    [15] 

Jo        if  no  A  solves  det{[I  +  MA]}  =  0 
[  else     (minA{cf(A)    |  det {[I  +  M A]}  =  O})'1 
While  this  expression  is  rather  unwieldy,  a  more  useful  alternative  formulation  is  given 

by 

n(Ma)  =  sup p(MaUd)  <  inf  aiDMaD-1).  (3.22) 

The  right-hand  side  of  Equation  3.22  has  been  proven  convex  (with  respect  to  scaling 

matrix  eD)    [28,  29]  and  the  MPDA  formulation  for  similarity  scaling  shows  that  the 

inequality  becomes  an  equality  at  the  "inf"  for  stationary  points  of  ^(DMaD'1)   [17]. 

Chapter  4  reviews  conditions  required  for  the  optimal  solution  to  Equation  3.22. 

However,  an  interesting  suboptimal  solution  was  proposed  by  Safonov  involving  the 

Frobenius  norm  minimization  of  DMaD~x     [14].     Since  the  Frobenius  norm  of  a 

matrix  provides  an  upper  bound  on  the  2-norm  (and  hence  the  maximum  singular 

value)  of  the  matrix,  reducing  the  Frobenius  norm  of  DMaD~l  generally  serves  to 

reduce  ^{DMaD'1)  as  well.  An  iterative  algorithm  introduced  by  Osborne  performs 

such  a  reduction  by  choosing  elements  of  D  that  equalize  the  row  and  column  sums  of 

DMaD~l   [30].  While  this  procedure  produces  a  conservative  upper  bound  for  //(M0), 

it  is  numerically  inexpensive  since  no  singular  value  decompositions  are  required. 
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3.4     Block  Structured  Uncertainties  with  Similarity  Scaling 

As  discussed  in  Section  2.4,  the  family  of  block  structured  uncertainties  represents 
the  most  general  class  since  it  is  a  superset  of  both  the  unstructured  and  element- 
by-element  structured  uncertainty  classes.  Let  A;,  represent  a  block  diagonal  matrix 
containing  both  scalar  and  nonscalar  uncertainty  elements.  As  in  Section  3.3,  this 
can  be  decomposed  giving 

Afe  =  PbUb 

where  Pb  contains  the  magnitude  bounds  and  Ub  the  phase  information  of  At.  An 
important  difference  between  this  decomposition  and  that  of  Section  3.3  is  the  re- 
quirement for  Pb  and  Ub  to  reflect  the  block  structure  of  A&.  For  this  configuration, 
the  diagonal  scaling  matrix  Db  must  also  reflect  the  block  structure  requiring  repeated 
elements  of  Dj,  as  necessary  to  correspond  to  the  size  of  the  individual  blocks  in  At. 
With  the  exception  of  maintaining  the  block  structure,  this  formulation  is  identical 
to  that  for  element-by-element  uncertainties  and,  in  fact,  the  block  form  of  Equation 
3.22  simply  becomes 

ft(Mb)  =  sup p(MbUb)  <  inf  a{DhMhDbl )  (3.23) 

Ub  Db 

where  Mb  =  MPb. 

3.5     Diagonalizing  Transformations 

As  mentioned  earlier,  it  is  always  possible  to  reconfigure  any  structured  A  into 
diagonalized  form.  First  consider  uncertainty  class  Ds  described  in  Section  2.4  where 
A   6  CnXn  contains  n2  complex  elements  St  to  6ni.     Next  define  A^   (E   Cn  Xn     = 
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diag[Sx, . . . ,  Sn2].  It  is  always  possible  to  relate  Ad  to  A  using  transformation  matrices 
Ex  G  $tnXn2  and  E2  €  9£n2xn  containing  only  l'j  and  0's  such  that  E^dE2  =  A.  From 
this  relationship,  the  necessary  and  sufficient  stability  conditions  of  Equation  3.1  may 
be  written  as 

supp(MExAdE2)  <  1. 

At  this  point,  some  means  of  recovering  the  simplifying  properties  of  diagonal  matrices 
required  by  Equation  3.20  is  needed.   The  following  theorem  relates  the  eigenvalues 
(and  hence  the  spectral  radius)  of  matrix  products  (AB)  and  (BA). 
Theorem  3.3:  Given  A  €  CnXm  and  B  G  CmXn, 

Xi(AB)  =  Xi(BA)     i  =  l,...,n. 

Proof:  Define  A4-  and  \j  as  the  eigenvalues  of  (AB)  and  (BA)  respectively.  Then 

(AB)Wi  =  XiWi     i  =  l,...,n 
(BA)^  =  A^     i  =  l,...,m 

where  Wi  G  Cnxl  and  Yj  G  CmXl  form  the  corresponding  eigenvectors.  Multiply  the 
first  equation  by  B  on  the  left  giving 

(BAB)Wi  =  XiBWi. 

Vector  BWi  G  CmXl  must  be  an  eigenvector  corresponding  to  the  ith  eigenvalue  of 
(-BA)  as  well  as  the  ith  eigenvalue  of  (AB).  Since  the  eigenvalues  of  a  matrix  are 
unique,  Xi(AB)  =  Xi(BA)  for  i  =  1, . . .  ,n.  ■ 
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Using  this  result,  replace  A  by  MExPd  and  B  by  E2  so  the  spectral  radius  equa- 
tions become 

p(MA)  =  p(E2ME1Ad). 

Next,  define  Ma  =  E2MEi  and  the  necessary  and  sufficient  stability  condition  of 
Equation  3.1  may  be  written 

supp(MaAd)  <  1. 

Starting  at  this  point,  the  analysis  of  Section  3.3  can  then  continue  as  before. 

While  this  diagonalizing  transformation  is  straightforward,  it  has  the  unfortunate 
consequence  of  increasing  the  system  size  from  n  x  n  to  as  much  as  n2  x  n2  with  a 
resultant  increase  in  computational  requirements  to  compute  the  maximum  singular 
values.  Also,  the  number  of  optimization  variables  required  to  infimize  Equation  3.22 
can  increase  to  n2  —  1  versus  2(ra  —  1)  for  nonsimilarity  scaling.  (As  with  matrices  L 
and  R  of  Equation  3.17,  a  single  element  of  D  may  be  arbitrarily  fixed  allowing  the 
elimination  of  one  optimization  variable).  Even  with  these  disadvantages,  the  ability 
to  manage  block  structured  uncertainties  makes  similarity  scaling  an  invaluable  tool 
in  the  analysis  of  uncertain  systems. 


CHAPTER  4 
THE  MAJOR  PRINCIPAL  DIRECTION  ALIGNMENT  PROPERTY 

4.1     Application  to  Similarity  Scaling 

The  use  of  singular  value  techniques  in  the  form  of  nonsimilarity  and  similarity 
scaling  has  the  disadvantage  of  initially  reducing  the  necessary  and  sufficient  stabil- 
ity conditions  of  Equation  3.1  to  simply  sufficient  conditions.  For  convenience,  the 
relevant  equations  are  repeated. 

supp(MA)     <    Ma{R-xML-l)a{LPR)     Nonsimilarity  Scaling         (4.1) 
snpp(MA)    <    mia(DMaD-1)     Similarity  Scaling.  (4.2) 

Eliminating  conservatism  in  these  singular  value  formulations  requires  establishing 
conditions  for  which  the  inequalities  hold  with  equality.  This  section  reveals  that 
achieving  equality  in  Equations  4.1  and  4.2  is  simply  a  matter  of  invoking  the  MPDA 
property  [17].  Although  MPDA  conditions  have  been  shown  to  exist  for  both  non- 
similarity  and  similarity  formulations,  this  discussion  will  focus  only  the  similarity 
scaling  techniques. 

Definition  4.1:  The  major  right  principal  direction,  Yj,  and  the  major  left  prin- 
cipal direction,  X1,  of  matrix  A  are  the  right  and  left  singular  vectors  of  matrix  A 
corresponding  to  a(A). 
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Theorem  4.1:  For  any  matrix  A  €  CnXn 

p(A)  =  a(A) 

if  and  only  if  X\  and  Y\  are  aligned  within  a  scaling  factor  eje  such  that 

Y1  =  e36X1.  (4.3) 

Proof:  For  notational  simplicity,  denote  a  (A)  by  a  and  p(A)  by  p.  Vectors  Xi  and 
Yi  from  the  singular  value  decomposition  must  be  unique  with  respect  to  one  another 
within  a  scaling  factor  e^.  Multiplying  both  sides  of  Equation  4.3  on  the  left  by  A 
gives 

AYX  =  e?6AXx. 

Also,  the  relationships  of  Section  3.1  show  that  AYX  =  ~&X\  so  this  can  be  rewritten 
as 

AXX  =  ae~seX1 

indicating  that  ae~jS  actually  corresponds  to  an  eigenvalue  of  A.  Since  eigenvalues 
of  a  matrix  are  always  magnitude  bounded  by  the  maximum  singular  value  of  the 
matrix,  sufficiency  of  the  theorem  is  established. 

Establishing  the  necessity  of  the  theorem  starts  with  the  assumption  that  a  —  p 
so  some  eigenvector  Z\  exists  where 

AZi  =  e*VZi  (4.4) 

and 

Z*1  AH  =  e'^uZf.  (4.5) 
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Multiplying  Equation  4.4  on  the  left  by  Equation  4.5  and  noting  the  AH  A  =  YY?Y 
provides  the  following  expression: 

ZjHAHAZ1       Z1HYY?YHZ1  _  __2 
ZfZx      "    Z?YY*ZX     ~  °  ' 

This  may  be  simplified  by  defining  a  new  vector  W\  =  YxH Z\  to  give 

WXHY?WX      _2  ,      , 

-w*wT  =  a  '  (4,6) 

Equation  4.6  can  only  be  satisfied  if  Wx  =  e3^ex  where  ex  is  the  first  standard  basis 
vector.  Therefore,  Zx  =  YWX  =  e^Yex  =  e^Yx.  Substituting  back  into  Equation  4.4 
then  gives 

Ae^Yt  =  We^Xi  =  e^e^Yx 

or 

Xx  =  e^Yx 

establishing  the  necessity  of  the  theorem.  ■ 

Corollary  4.1:  For  the  case  of  repeated  singular  values  of  A  6  CnXn  where  q  denotes 
the  multiplicity  of  ax, 

p(A)  =  ax(A),...,aq(A) 

if  and  only  if  the  subspaces  spanned  by  Xx , . . . ,  Xq  and  Yx,. . .  ,Yq  are  aligned  within 
a  scaling  factor  ej6 . 

Proof:  From  Equation  3.2  the  singular  vectors  of  matrix  A  with  a  repeated  maxi- 
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mum  singular  value  are  denned  as 

AHAY1  =  a\Yx  AAHXX  =  a\Xx 

\  and  : 

AHAYq  =  a\Yq  AAHXq  =  u\Xq 

Therefore,  linear  combinations  of  the  first  q  columns  of  X  or  Y  form  legitimate 
solutions  for  X\  and  Yt  respectively.  Combining  these  solutions  with  Theorem  4.1 
completes  the  proof.  ■ 

The  MPDA  property  therefore  provides  an  analytical  test  for  equality  in  Equation 
3.21  through  which  the  necessary  and  sufficient  stability  conditions  of  Equation  3.1 
apply.  This  allows  the  innmization  of  Equation  3.22  to  be  performed  while  providing 
a  direct  test  as  to  whether  the  structured  singular  value,  (//(Ma)),  has  been  achieved. 
It  would  be  desirable  if  this  innmization  was  somehow  guaranteed  to  invoke  MPDA 
and  hence  attain  /J,(Ma).  Fortunately,  such  a  guarantee  is  provided  by  the  following 
theorem. 

Theorem  4.2:  The  infimizing  solution  with  respect  to  diagonal  scaling  matrix  D 
of  stability  criterion 

inf  aiDMaD'1)  <  1  (4.7) 

is  both  necessary  and  sufficient  provided  that  J^  =  0  at  the  "inf." 

Proof:  The  sufficiency  of  the  theorem  was  shown  in  Equation  3.21.  Necessity  is 
established  as  long  as  MPDA  can  be  achieved.  For  simple  a(DMa_D_1),  the  optimal 
solution  of  Equation  4.7  occurs  when 

da2 


ddi 


(DMaD-1)  =  0 
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indicating  a  stationary  point.  As  mentioned  in  Section  3.3,  the  left-hand  side  of 
Equation  4.7  is  convex  with  respect  to  scaling  matrix  eD  so  any  local  minimum  is  the 
global  minimum.  Therefore  the  only  stationary  point  must  correspond  to  the  optimal 
solution.  By  direct  differentiation,  the  stationarity  condition  becomes 

ffcj2  2YH 

■jfc (Ma)  =  -±-  [M»EiMa  -  a2Ei]  Yi  =  0  (4.8) 

where  Ma  =  MaD_1,  and  E\  is  a  square  matrix  containing  a  "1"  in  the  nth  position 
and  "Os"  elsewhere.  (Note:  the  actual  differentiation  is  omitted  here  but  the  complete 
differentiation  of  a  similar  function  appears  in  Section  5.5). 

Examining  conditions  at  which  Equation  4.8  equals  zero  reveals  the  requirement 
that 

l*i<|  =  Iviil     i  =  l,...,n  (4.9) 

which  is  in  fact  the  MPDA  requirement  that  the  magnitudes  of  the  elements  of  Xt 
and  Y\  be  equal. 

Completion  of  the  proof  from  this  point  requires  that  some  Ud  exists  such  that 
the  phase  alignment  of  Xx  and  Yx  occurs.  The  existence  of  such  a  Ud  then  guarantees 
that  p(MaUd)  =  a(DMaD~lUd)  =  n{Ma).  From  the  singular  value  decomposition, 
the  principal  singular  vectors  are  related  as 

X?(DMaD~l)    =    aY*. 
Introducing  a  diagonal  unitary  matrix,  Ud,  allows  these  equations  to  be  rewritten  as 


37 


{DMaD~lUd)UdHY1    =    aX1 
X*{DMaD-lUd)    =    oY?Ud.  (4.10) 

From  Equation  4.10  it  is  apparent  that  the  diagonal  elements  of  Ud  may  be  used  to 
alter  the  phases  of  Xi  and  Y\  while  continuing  to  satisfy  the  MPDA  requirement  that 
the  magnitudes  of  the  corresponding  elements  of  Xi  and  Yx  be  equal.  The  ability  to 
achieve  independent  phase  alignment  of  Xx  and  Yx  by  Ud  completes  the  proof.  ■ 

The  result  of  this  theorem  is  that  optimization  of  Equation  4.7  has  the  effect 
of  directly  invoking  MPDA  at  the  "inf"  for  simple  a,  guaranteeing  the  necessary 
and  sufficient  stability  conditions.  Cases  with  nonsimple  a  (known  as  "cusps")  do 
not  necessarily  contain  a  stationary  point  at  the  "inf"  and  must  be  analyzed  using 
different  techniques.  A  more  complete  discussion  of  cusping  systems  appears  in 
Chapter  8. 

4.2     Block  Structured  Uncertainties 

By  following  the  MPDA  development  for  scalar  uncertainties,  a  similar  formu- 
lation results  for  general  block  structured  uncertainties.  The  main  difference  is  the 
requirement  of  maintaining  the  individual  block  structures  rather  than  considering 
each  element  of  A  individually.  As  discussed  in  Section  3.4,  this  block  structure  must 
be  reflected  in  the  scaling  matrices  as  well.  The  optimal  solution  is  then  characterized 
by  the  equalization  of  the  2-norms  in  the  corresponding  blocks  of  the  left  and  right 
principal  vectors  X%,  and  Y\  respectively.  The  same  guarantee  of  achieving  equality 
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at  the  optimal  solution  of  Equation  3.23  for  simple  a(DbMhD'l;1)  applies  as  in  the 
element-by-element  structured  uncertainty  case. 

4.3     Examples 


The  following  examples  illustrate  the  methods  described  in  the  previous 


sec- 


tions. 


Example  4.1:  Let  randomly  generated  matrices  M  £  C3x3  be  chosen  such  that 


M 


.063 +  .156J  -.322 +.480;  .585 +  .526; 
.726  -.514;  -.323  -.344;  .150  -  .469; 
.189 -.463;       .053  -  .577;    -.236  -  .056; 


Next  chose  a  random  element-by-element  structured  uncertainty  matrix  A  £  C3x3 
with  corresponding  P  £  9£3x3  written  as 


N  l*a  I  N 
M  N  |4| 
N   N   N 


< 


2.99  3.03  0.54 
1.65  1.87  3.41 
1.90    1.20    1.37 


Formulated  as  a  nonsimilarity  scaling  problem, 


fi(MA)  <  miaiR^ML^piLPR) 


the  diagonal  elements  of  the  optimal  R  and  L  scaling  matrices  are  found  to  be 


R 


The  maximum  singular  values  do  not  repeat  for  this  example  so 


n  =  i 

r/1  =  l 

r2  =  1.151 

,  L  = 

h  =  0.944 

r3  =  1.174 

?a  =  1-203 

fi(Ma)  =  a(R-1ML~1)a(LPR)  =  8.25. 

This  same  problem  may  be  recast  into  a  similarity  scaling  form  by  first  selecting 
transformation  matrices  E\  and  E2  as 
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£i  = 


1110  0  0  0  0  0 
0  0  0  1110  0  0 
0    0    0    0    0    0    111 


"10    0" 

0    1    0 

0    0    1 

1    0    0 

,  £2  = 

0    1    0 

0    0    1 

1    0    0 

0    1    0 

.001. 

Expanding  M  and  P  as  described  in  Section  3.3,  the  optimal  solution  of  Equation 

3.22  requires  scaling  matrix  D  with  diagonal  elements 

di  =  1 
d2  =  0.906 
4  =  0.406 
d4  =  0.719 
dh  =  0.688 
d6  =  0.988 
d7  =  0.783 
4  -  0.560 
d9  =  0.636 

As  before,  n(Ma)  =  ^(DMaD-1)  =  8.25.   The  absolute  values  of  the  left  and  right 
principal  directions  for  this  optimal  solution  are 


1*1 


1 

'  1 

1.119 

1.119 

0.444 

0.444 

0.719 

0.719 

0.850 

=  \Yt\  = 

0.850 

1.080 

1.080 

0.783 

0.783 

0.691 

0.691 

0.695  . 

0.695 

As  required  by  the  MPDA  theory,  [Xi|  =  \Yi\  at  the  "inf"  guaranteeing  that  some 
Ud  exists  such  that  p{MaUd)  =  a(DMai)-1)  =  n{Ma). 

Example  4.2:  Using  the  M  matrix  from  Example  4.1,  define  a  block  structured 
uncertainty  matrix  as 
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A  = 


Si   s2   s3 

85       A6 


where  81  through  £5  are  scalar  uncertainties  such  that  |£,-|  <  p,-  and  Ab  is  a  2  x  2  block 
with  a(Ab)  <  pb.  Block  Ab  may  be  represented  as 

A6  =  PbUb 

where  Pb  is  a  2  x  2  matrix  with  a(A)  =  J?j  and  £4  is  a  2  x  2  unitary  matrix  containing 
the  phase  information  of  Ab.  We  may  now  define  a  matrix  PA  containing  randomly 
generated  bounds  on  the  8{  and  Ab  elements  as 


Pa 


2.99    3.03    0.54 

1.65 

1.90  Ph 


with  a(Pb)  =  1.87.  It  is  easily  shown  that  the  worst  case  structure  for  Pb  in  terms  of 
/i(M;,)  occurs  when 


Ph  = 


1.87      0 
0      1.87 


allowing  a  P  matrix  for  this  case  to  be  written  as 

P  = 


2.99    3.03    0.54 
1.65    1.87      0 
1.90      0      1.87 


The  2x2  block  structure  constrains  the  last  two  terms  of  scaling  matrix  Dh  to  be 
identical  thus  eliminating  one  degree  of  freedom  from  the  optimization.  Also,  the  two 
zero  elements  of  P  allow  the  diagonalized  matrix  Pd  (and  hence  Mb)  to  be  of  size 
7x7  rather  than  9x9. 


Choose 
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Ei  = 


1110  0  0  0 
0  0  0  10  10 
0    0    0    0    10    1 


"  1    0    0  " 

0    1    0 

0    0    1 

,  E2  = 

1    0    0 

1    0    0 

0    1    0 

0    0    1 

corresponding  to  Pd  =  diag[2M,  3.03,  0.54,  1.65,  1.9,  1.87,  1.87]  so  that  P  = 
-EiPd£2.  Expanded  matrix  Mb  may  then  be  formed  as  Mb  =  E2ME1Pd.  Following 
the  infimization  of  a{D%MhD^%  the  optimizing  Db  and  the  absolute  values  of  the 
elements  of  Xx  and  Y^  are  found  to  be 


A 


rfi  =  1 
d2  =  0.82 
d3  =  0.41 
d4  =  0.63 
d5  =  0.87 
<k  =  0.69 
d7  =  0.69 


1-^x1  = 


0.427 
0.527 
0.187 
0.269 
0.373 
0.442 
0.313 


1511  = 


0.427 
0.527 
0.187 
0.269 
0.373 
0.279 
0.465 


with  fi(Mb)  =  6.5.  Note  that  the  first  five  terms  of  \XX\  and  \y\\  are  equal  as  in  the 
case  of  scalar  uncertainty  blocks  under  MPDA  conditions.  However,  the  presence  of 
the  2  x  2  block  causes  the  last  two  terms  of  \Xj  |  and  \YX\  to  be  related  in  the  following 
manner: 

l^eil2  +  |x71|2  =  |y61|2  +  \yn\2  =  0.2937. 

As  described  in  Section  4.2,  this  equality  between  norms  of  the  corresponding  blocks  of 
Xj  and  Yi  occurs  whenever  MPDA  is  established  with  block  structured  uncertainties 
guaranteeing  the  existence  of  some  optimal  Ub  such  that  p(MbUb)  =  ^(DhMht)^1)  = 
f*(Mb)    [31]. 

The  material  presented  so  far  provides  a  motivation  for  computing  the  structured 
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singular  value  using  singular  value  scaling  techniques.  However,  even  with  the  ex- 
amples, very  little  insight  can  be  gained  into  the  computational  issues  involved  with 
determining  ft.  In  fact,  although  the  structured  singular  value  concept  is  currently 
being  applied  to  the  design  of  large  engineering  problems  [23],  a  number  of  difficul- 
ties continue  to  hinder  its  application.  The  next  chapter  shows  that  knowledge  of  the 
uncertainty  structure  can  lead  to  a  reduction  in  the  number  of  optimization  variables 
required  to  compute  \i  using  similarity  scaling. 


CHAPTER  5 

REDUCTION  IN  THE  NUMBER  OF  OPTIMIZATION  VARIABLES  REQUIRED 

TO  COMPUTE  THE  STRUCTURED  SINGULAR  VALUE 

5.1     Reduction  of  Optimization  Variables 

As  mentioned  earlier,  similarity  scaling  techniques  have  the  capability  of  handling 
general  block  structured  uncertainties.  While  this  advantage  justifies  the  use  of  simi- 
larity scaling  for  systems  with  these  uncertainty  structures,  the  requirement  that  the 
uncertainty  blocks  be  diagonalized  results  in  a  possibly  substantial  increase  in  overall 
system  size.  For  scalar  blocks  this  expansion  transforms  annxn  system  to  an  n2  x  n2 
system.  Not  only  does  this  extended  system  size  increase  the  floating  point  operation 
(flop)  count  for  each  singular  value  decomposition,  but,  as  discussed  in  Section  3.5, 
it  also  increases  the  number  of  optimization  variables  from  2(n  —  1)  for  nonsimilarity 
scaling  to  n2  —  1  for  similarity  scaling.  Since  M,  A  and  hence  (i(Ma)  are  frequency 
dependent  functions,  a  frequency  sweep  of  fj,  must  be  performed  to  determine  the 
worst  case  condition.  Depending  on  the  operating  frequency  range  of  the  control 
system  under  investigation,  this  frequency  sweep  may  involve  hundreds  or  thousands 
of  individual  points  at  which  /J,  must  be  evaluated.  Each  of  these  points  therefore 
represents  a  substantial  cost  in  CPU  time. 

Since  both  nonsimilarity  and  similarity  scaling  methods  produce  identical  results 
for  element-by-element  structured  uncertainties,  it  would  seem  reasonable  to  con- 
clude that  some  of  the  n2  —  1  degrees  of  freedom  required  for  similarity  scaling  are 
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actually  redundant.  Removing  these  unnecessary  optimization  variables,  if  possible, 
should  provide  for  some  corresponding  reduction  in  flops  required  to  compute  fi  using 
similarity  scaling.  The  following  theorem  and  subsequent  proof  show  that,  in  fact, 
no  more  that  2(ra  -  1)  optimization  variables  are  required  to  compute  fj,  for  both 
similarity  and  nonsimilarity  scaling. 

Theorem  5.1:  Given  M  e  CnXn  and  A  <=  CnXn,  where  A  is  composed  of  1  x  1 
blocks,  the  solution  for  fi(Ma)  in  the  optimization  problem 

KMa)<Ma(DMaD-1)  (5.1) 

can  be  found  with  no  more  than  2(ra  -  1)  free  variables  as  long  as  a  stationary  point 
occurs  at  the  "inf." 

Proof:  To  reduce  notational  complexity  this  proof  will  be  developed  explicitly  for 
the  2  x  2  case  and  the  general  n  x  n  case  will  follow  as  a  simple  extension.  Let 


M  = 


mn     m12 
m2i     m22 


\Si\   |4| 


< 


Pn    Pu 

P21      P22 


8X    62 


and  Pd  =  diag\pn,p12,jy21,p22].  For  this  Pd  choose 

E1  = 


'10  1 

110    0 
0    0    11 

,    E2  = 

0  1 

1  0 

.01. 

Applying  scaling  matrix  D  =  diag{dx,  d2,  d3,  d4]  to  Ma  gives 


DMaD~l  = 


mnPu        %rnnPl2    Jm12p21     ^m12p22 


^m2ipn  m2lPl2  f3m22p21  fm22p22 
%™-\\V\\  J"*iiPi2  rn12p2Y  fml2p22 
^mnpn     %rn2XpX2    Jm22p21     m22p22 


(5.2) 
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Also,  from  the  singular  value  decomposition, 


DMaD'1  =  XVYH 


(5.3) 


where  the  orthogonal  matrices  X  and  Y  represent  the  left  and  right  singular  vectors 
of  (DMaD-x)  respectively  and  E  is  a  positive  diagonal  matrix  containing  the  singular 
values  of  (DMaD'1)  in  decreasing  order  of  magnitude. 

For  the  2x2  case,  matrix  M  must  have  rank<  2.  Multiplying  M  on  the  left  by 
E2  and  on  the  right  by  ExPd  to  form  Ma  therefore  requires  that  rank(Ma)<  2  forcing 
03  =  o"4  =  0.  Expansion  of  Equation  5.3  into  its  individual  elements  gives 


DMaD-x  = 


ZiiI7ii0-i  +  x12y12a2  ■  ■  •  xny41a-L  +  x12y42a2 

X2iyno-i  +  x22y12a2  ■  ■  ■  x21y41ai  +  x22y42a2 

£3i]/iicri  +  x32y12a2  ■  ■  ■  x31]/41cri  +  x32y42a2 

^4iyn0i  +  x42y12a2  •  •  •  x41y41a1  +  xi2y42a2 


(5.4) 


Examination  of  Equation  5.2  shows  that  the  ratio  of  the  first  and  third  rows 
is  a  scalar  constant  ax  =  o?i/d3.  Since  Equation  5.2  and  Equation  5.4  both  equal 
DMaD~x,  the  ratio  of  the  first  and  third  rows  of  Equation  5.4  must  also  be  ax. 
Equating  the  first  row  of  Equation  5.4  to  cxx  times  the  third  row  requires  that 


2/ii  = 

F21  = 

?31  = 

#41  = 


__      ^2(^12  -  <^\X32) 
V\2       7  v 

__   0-2(^12  -  0:1X32) 

2/22        "1  v 

<j\[xxi  -  axxzi) 

_    a2(x12  -  Q1X32) 

2/32     7  r 

cr^xn  -  0:1X31) 

_     gj(^12  -  Ql^32) 

2/42       7  \  ■ 

<Tl(Xn  -  0:1X31) 
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This  set  of  equations  implies  a  relationship  between  Yx  and  Y2  of  the  form  Yx  =  kY2 
where  k  =  ^j|™^j.  However  a  linear  dependence  between  Y1  and  Y2  is  impossible 
because  the  columns  and  rows  of  Y  are  orthogonal  to  one  another.  This  inconsistency 
is  satisfied  only  for  xu  =  axx31  and  x12  =  a\xZ2-  Relationships  between  the  remaining 
elements  of  Xx  and  Yt  may  be  developed  through  a  similar  process.  The  second  and 
fourth  rows  of  5.4  have  a  ratio  a2  =  d2/d4  requiring  that  x21  =  a2x41  and  x22  =  a2x42. 
Finally,  the  ratio  of  the  first  and  second  columns  of  Equation  5.4  is  a3  =  #£ii  and 

^  dlPl2 

the  ratio  of  the  third  and  fourth  columns  of  Equation  5.4  is  a4  =  ^221.  These 
relationships  between  the  columns  of  Equation  5.4  require  that  yn  =  a3y21  and 
2/3i  =  <*42/4i-  Summarizing  these  results  gives 


Xx  = 


~ 

- 

xn 

Vn 

X21 

,    Yx  = 

2/ll/«3 

Xn/ati 

2/31 

x2\/a2 

2/31/0:4 

(5.5) 


so  that  only  two  independent  terms  (xn,  x21  in  Xu  and  yn,  y3l  in  Yt)  appear  in 
each  vector.  From  Section  4.1,  equality  between  the  left-hand  and  right-hand  sides 
of  Equation  5.1  occurs  only  when  MPDA  conditions  are  established.  This  in  turn 
requires  pair-wise  equality  between  the  elements  of  Xx  and  Y\  or 


kill  =  \yn\    i  =   l,...,n2 


(5.6) 


for  similarity  scaling  with  simple  ^DMaD'1).  Applying  the  requirements  of 
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Equation  5.6  to  Equation  5.5  gives 


knl 

= 

toil 

k2i| 

= 

ls/xi| 

«3 

.  knl 

«3 

kn| 

<*1 

= 

12/31 1 

kail 



W 

knl 

toil 

a2  o4         Oi04       «2«3 

From  this  scries  of  relationships  the  following  equality  must  hold  under  MPDA  con- 
ditions: 

0:10:4 

o^  =  L  (5-7) 

Replacing  the  a*a  in  Equation  5.7  with  the  corresponding  d»a  and  Pij,s  gives 

djPuP2i  4 

which  shows  that  J4  can  be  expressed  as  a  simple  function  of  du  d2,  d3  and  the  Pij,s. 
The  proof  for  the  2  x  2  case  is  completed  by  noting  that  one  of  the  dt  can  be  always  be 
scaled  to  "1"  without  loss  of  generality  so  that  the  total  number  of  free  optimization 
variables  is  2  =  2(n  —  1). 

Having  completed  the  proof  for  the  2  x  2  case,  the  proof  for  the  general  n  x  n 
case  follows  in  the  same  manner.  For  M  G  CnXn  the  rank  of  expanded  matrix  Ma 
must  be  less  than  or  equal  to  n  requiring  that  <rn+u  .  ..,an2=0.  Continuing  with  the 
process  of  equating  rows  and  columns  of  Equations  5.2  and  5.4  produces  relationships 
between  Xx  and  Y1  similar  to  those  of  Equation  5.5  when  MPDA  is  established.  The 
number  of  independent  vector  elements  in  X1  and  Yx  becomes  n  for  the  general  case 
and  (n-1)2  expressions  in  the  form  of  Equation  5.7  describe  the  relationships  between 
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the  2n(n  -  1)  aitg.  From  these  expressions  a  general  equation  relating  the  dependent 
and  independent  values  of  optimal  scaling  matrix  D  can  be  written  as 


dlP(l,(i-n))P2,l 

i    =    n  +  2, . . .  ,2n; 

j2  d2j-2nd2n+lPnP3,(j-2n) 

3  dlP(l,U-2n))P3,l 

j     =    2n  +  2, . . . ,  3n; 


(5.8) 


4 


dk-(n-l)nd(n-l)n+lPnP(n,k-(n-l)n) 
d\P(l,k-{n-\)n)Pn,\ 


k    =    (n-l)n  +  2, 


,n 


Subtracting  the  (n  -  l)2  expressions  in  the  form  of  Equation  5.7  from  the  n2  -  1 
variables  required  for  the  unreduced  scaling  structure  gives  the  desired  result  of 


n 


l-(n-l)2  =  2(n-l) 


optimization  variables. 


The  complexity  of  the  indices  in  Equation  5.8  unfortunately  obscures  the  simplicity 
of  its  application.  The  following  example  will  hopefully  clarify  these  results. 

Example  5.1:  Let  randomly  generated  matrices  M  G  C3*3  and  P  G  3£3x3  be  chosen 

as  in  Example  4.1  such  that 

".063 +  .156.7    -.322  +  .480j        .585  +  .526; 

M=      .726 -.514;    -.323-.344j       .150  -  .469; 

_  .189  -  .463;       .053  -  .577;    -.236  -  .056; 


and 


P 


P\\  Pl2  Pl3 
P21  P22  P23 
P31       P32      P33 


2.99  3.03  0.54 
1.65  1.87  3.41 
1.90    1.20    1.37 


49 


Then  expanding  Ma  =  E2ME1Pd  gives  p,{Ma)  =  8.25  with  optimal  scaling  matrix 


D  = 


dx  =  1 
d2  =  0.906 
d3  =  0.406 
d4  =  0.719 

d5  =  0.688 

d6  =  0.988 
d7  =  0.783 
4  =  0.560 

dQ  =  0.636 


_      ld%diPiiP22 

d'fpi2P21 


<^?Pl3P21 

<*2<*7PllP32 
^lPl2P31 

^3^?PllP33 
^lPl3P31 


Of  course  this  is  identical  to  scaling  matrix  D  from  Example  4.1  found  by  optimizing 
over  n2  -  1  variables.  In  this  case,  the  added  complexity  of  computing  the  dependent 
elements  of  D  is  more  than  offset  by  the  savings  in  computations  resulting  from  the 
reduced  number  of  optimization  variables. 

Tables  5.1  through  5.3  summarize  the  savings  in  floating  point  operations  (flops) 
due  to  the  scaling  method  of  this  section.  The  "Size"  columns  represent  the  dimension 
of  expanded  matrix  Ma  while  the  number  of  optimization  variables  (n2-l,  or  2(n-l)) 
column  distinguishes  between  the  old  and  new  scaling  methods  respectively.  Both 
optimization  methods  were  implemented  using  a  MATLAB™  based  quasi-Newton 
method  with  analytic  derivatives  [32,  33].  Also,  the  starting  point  D  =  I  was  chosen 
for  both  methods  to  make  the  comparison  as  equitable  as  possible.  This  initial  choice 
for  D  is  probably  as  good  as  any  other  for  randomly  generated  M  and  P  matrices. 
However,  stability  analysis  of  an  actual  system  using  a  frequency  sweep  could  use  the 
optimum  D  from  a  previous  frequency  point  as  an  improved  initial  point. 

For  the  4x4  systems,  the  average  flop  count  required  to  compute  (i  using  n2  -  1 
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Table  5.1 

:  FLO 

P  Count  ( 

comparison  for 

1x4  Systems 

System 

Size 

Optimization  Variables 

%  Difference 

n2-l 

2(ra-l) 

1 

4 

85,847 

70,669 

-17.7 

2 

4 

84,463 

55,117 

-34.7 

3 

4 

102,588 

95,816 

-6.6 

4 

4 

83,236 

74,326 

-10.7 

5 

4 

83,239 

62,915 

-24.4 

6 

4 

67,001 

57,630 

-14.0 

7 

4 

75,247 

61,280 

-18.6 

8 

4 

62,237 

73,254 

+17.7 

9 

4 

78,078 

43,474 

-43.6 

10 

4 

61,244 

58,959 

-3.7 

11 

4 

77,018 

34,567 

-55.1 

12 

4 

51,155 

48,680 

-4.8 

13 

4 

88,406 

73,190 

-17.2 

14 

4 

53,459 

46,563 

-14.7 

15 

4 

57,321 

46,417 

-19.0 

16 

4 

59,942 

49,500 

-17.4 

17 

4 

47,572 

33,124 

-30.3 

18 

4 

56,281 

51,606 

-8.3 

19 

4 

76,508 

61,132 

-20.1 

20 

4 

89,680 

54,088 

-39.7 

21 

4 

63,493 

68,793 

+8.3 

22 

4 

75,972 

58,982 

-22.4 

23 

4 

74,279 

101,329 

+36.4 

24 

4 

77,029 

53,882 

-30.04 

25 

4 

76,402 

39,261 

-48.6 

Average: 

75,637 

57,955 

-23.4 
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Table  5.2:  FLO 


System      Size 


Count  Comparison  for  9  x  9  Systems 


Optimization  Variables 


n2  -  2(ra-l) 


%  Difference 


1 

9 

1,178,572 

1,240,971 

5.3 

2 

9 

1,188,389 

600,863 

-49.4 

3 

9 

1,298,812 

991,019 

-23.7 

4 

9 

1,079,390 

810,888 

-24.9 

5 

9 

1,146,836 

785,305 

-31.5 

6 

9 

1,924,664 

1,315,678 

-31.6 

7 

9 

2,320,320 

1,172,074 

-49.5 

8 

9 

1,596,230 

1,332,371 

-16.5 

9 

9 

1,356,005 

1,182,523 

-12.8 

10 

9 

1,591,708 

1,294,051 

-18.7 

11 

9 

772  825 

606,518 

-21.5 

12 

9 

784,541 

922,383 

+17.5 

13 

9 

897,372 

798,991 

-11.0 

14 

9 

1,093,402 

1,090,522 

-5.8 

15 

9 

906,052 

501,610 

-44.6 

16 

9 

824,228 

847,159 

+2.8 

17 

9 

820,406 

731,396 

-10.8 

18 

9 

1,045,708 

865,101 

-17.3 

19 

9 

897,115 

605,850 

-32.5 

20 

9 

1,538,927 

982,515 

-36.2 

21 

9 

1,366,940 

938,130 

-31.5 

22 

9 

932,012 

722,522 

-22.5 

23 

9 

1,190,635 

581,242 

-51.2 

24 

9 

943,404 

990,477 

+4.9 

25 

9 

1,270,973 

831,299 

-17.5 

Average: 


1,167,769 


909,658 


-22.0 
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Table  5.3: 

FLOP  Count  Comparison  for  16  x  16  Systems 

System 

Size 

Optimization  Variables 

%  Difference 

n2-l 

2(n-l) 

1 

16 

8,124,858 

5,633,607 

-30.7 

2 

16 

9,016,172 

4,483,259 

-50.3 

3 

16 

10,886,765 

9,511,034 

-12.6 

4 

16 

8,247,534 

6,271,272 

-24.0 

5 

16 

8,519,069 

5,544,709 

-34.9 

6 

16 

9,450,057 

5,919,971 

-37.4 

7 

16 

10,789,909 

9,097,814 

-15.7 

8 

16 

9,574,741 

8,900,847 

-7.0 

9 

16 

10,896,836 

8,523,835 

-21.8 

10 

16 

11,425,524 

7,454,489 

-34.8 

11 

16 

8,525  125 

6,322,816 

-25.8 

12 

16 

8,514,427 

5,807,702 

-31.8 

13 

16 

9,030,860 

8,373,217 

-7.3 

14 

16 

7,637,353 

5,854,878 

-23.3 

15 

16 

10,589,478 

5,854,878 

-43.9 

16 

16 

10,629,956 

6,544,465 

-38.4 

17 

16 

9,661,623 

8,417,042 

-12.9 

18 

16 

9,175,808 

5,742,484 

-37.4 

19 

16 

7,354,233 

6,232,710 

-15.3 

20 

16 

7,611,841 

5,889,917 

-22.6 

21 

16 

10,710,469 

6,438,349 

-39.9 

22 

16 

11,214,534 

9,057,426 

-19.2 

23 

16 

7,954,635 

6,603,402 

-17.0 

24 

16 

9,350,555 

6,028,401 

-35.5 

25 

16 

8,541,532 

7,901,464 

-7.5 

Average: 

8,996,696 

6,896,279 

-23.3 
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variables  is  75,637  versus  57,955  using  2{n  -  1)  variables:   a  savings  of  23.4%.    As 

system  size  increases,  the  flop  count  savings  do  not  change  appreciably  with  an  average 
savings  of  22.0%  for  the  9  x  9  systems  and  23.3%  for  the  16  x  16  systems.  Since 
li  must  be  determined  over  a  wide  frequency  range  for  each  transfer  function,  the 
savings  offered  by  the  scaling  method  of  Theorem  5.1  promise  substantial  reductions 
in  computational  effort  for  the  design  of  robust,  multivariable  control  systems.  Note 
that  in  some  of  the  examples  shown  the  reduced  scaling  method  actually  requires 
more  flops  than  the  full  optimization  over  n2  -  1  variables.  This  occasional  anamoly 
is  probably  due  to  the  initial  guess  of  D  =  I.  For  the  unreduced  case,  all  diagonal 
elements  of  D  may  be  set  to  1,  while  the  reduced  structure  only  allows  the  free 
elements  to  be  set  to  1.  For  a  few  of  the  examples,  having  all  diagonal  elements  of 
D  =  1  apparently  provides  the  unreduced  method  with  a  superior  initial  guess  for 
the  optimization. 

Reducing  the  number  of  optimization  variables  offers  more  than  the  savings  in 
flop  counts.  Depending  on  the  minimization  algorithm  chosen,  a  possibly  substantial 
savings  in  memory  requirements  results  as  well.  For  example,  the  popular  quasi- 
Newton  (such  as  that  used  in  this  work)  routines  require  storage  on  the  order  of 
N2  where  N  represents  the  number  of  optimization  variables  [34].  Along  with  the 
memory  requirements  comes  additional  bookkeeping  to  keep  track  of  the  gradient 
information  for  the  N  variables.  Therefore  the  reduction  from  n2  -  1  to  2(n  -  1)  free 
variables  offers  computational  savings  that  may  not  show  up  directly  in  flop  count 
comparisons. 
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5.2     Effects  of  Zero  Elements  in  P 


Example  5.1  and  the  systems  of  Table  5.1  represent  the  maximum  reduction  of 
free  variables  from  n2  -  1  to  2(n  -  1).  If  elements  of  A  are  replaced  with  zeros, 
this  efficiency  gap  begins  to  diminish  although  the  new  method  is  always  at  least  as 
efficient  in  terms  of  optimization  variables  as  previous  similarity  scaling  methods.  In 
fact,  the  efficiency  advantage  holds  until  the  number  of  zeros  equals  n(n  -  2)  +  1  so, 
for  example,  a  10  x  10  system  with  up  to  80  zero  elements  in  A  will  still  maintain  an 
advantage  in  terms  of  free  variables.  Of  course  the  relationships  shown  in  Equation 
5.8  must  be  altered  to  account  for  the  reduction  in  the  dils.  This  is  demonstrated  in 
Example  5.2. 

Example  5.2:  Starting  with  the  same  random  M  and  P  matrices  from  Example 
4.1,  set  p22  and  p33  to  zero.  Then  the  size  of  expanded  matrix  Ma  need  only  increase 
to  7  x  7  instead  of  9  x  9  since  the  zero  elements  of  P  can  be  eliminated.  The  optimal 
D  for  this  example  is  then  found  to  be 


<k  =  1 

d2  =  0.852 
da  =  0.387 
d4  =  0.746 

d5  =  0.976 

d6  =  0.868 

d7  =  0.584 


D  = 


d'jPl3P21 


<%.<%P\\PZ2 


with  fj,(Ma)  =  6.64.  Note  that  while  the  expression  for  d5  corresponds  to  that  of  d6 
in  Example  5.1,  the  expression  for  d7  changes  slightly  to  account  for  the  reduction  in 
the  number  of  d{,s.  This  example  shows  the  reduction  from  n2  -  1  =  8  to  n2  -  3  =  6 
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for  previous  methods  of  similarity  scaling.  The  new  method,  however,  experiences  no 
additional  reduction  and  remains  at  2(n  —  1)  =  4 

5.3     Application  to  Block  Structured  Uncertainties 

As  discussed  in  Section  2.4,  block  structured  uncertainties  represent  an  impor- 
tant class  of  uncertainties  that,  in  general,  cannot  be  addressed  using  nonsimilarity 
scaling  techniques.  Therefore,  one  of  the  more  important  advantages  of  similarity 
scaling  is  the  ease  with  which  general  block  structured  uncertainties  may  be  handled. 
Furthermore,  by  examining  the  relationships  between  the  elements  of  Xx  and  Y\  as 
in  Equation  5.5,  it  is  also  possible  to  reduce  the  number  of  optimization  variables  for 
blocks  of  arbitrary  dimension.  Theorem  5.2  follows  directly  from  Theorem  5.1  and 
the  MPDA  relationships  for  block  structured  uncertainties  discussed  in  Section  4.2. 

Theorem  5.2:  Given  M  e  CnXn  and  A  €  CnXn,  where  A  is  composed  of  blocks  of 
dimension  less  than  n,  the  solution  for  n(Mb)  in  the  optimization  problem 

fi(Mb)<Ma(DbMaD^)  (5.9) 


with  simple  a  can  be  found  with  no  more  than  2(n  —  1)  free  variables. 

Proof:  The  proof  is  straightforward  given  Theorem  5.1.  Recall  from  Section  3.4 
that  for  block  structured  uncertainties,  scaling  matrix  Db  contains  block-diagonal 
elements  corresponding  to  the  structure  of  P.  For  P  composed  solely  of  1  x  1  blocks, 
Theorem  5.1  reveals  that  no  more  than  2(n  -  1)  optimization  variables  are  required 
to  compute  fi(Mb).  At  the  other  extreme,  for  P  composed  of  a  single,  diagonal,  block 
structured  uncertainty  of  size  n,  no  optimization  variables  would  be  required  since 
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ems 


the  system  effectively  simplifies  to  an  unstructured  uncertainty  problem.  For  syst 
falling  between  these  limiting  cases,  singular  vector  relationships  similar  to  those  of 
Equation  5.5  still  apply  so  the  number  of  optimization  variables  required  to  solve  any 
block  structured  problem  cannot  exceed  2(n  -  1).  ■ 

Example  5.3  illustrates  the  application  of  the  optimization  variable  reduction 
method  to  block  structured  uncertainties. 

Example  5.3:  Using  the  same  M  and  A  matrices  from  Example  4.2,  form  expanded 
matrix  Mb  followed  by  DMbD~\  The  rank  of  DMaD~*  must  be  less  than  or  equal 
to  3  so  comparing  the  terms  of  DMaD~l  in  the  form  of  Equation  5.2  with  that  in  the 
form  of  Equation  5.4  requires  the  structure  of  X1  and  Y±  to  be 


x1 


Xn 


^21 


^31 


Zll 


XU 


di 


X^i 


231 


yn 


"npnd 


•2 


Yx 


tfllPllrf3 


2/41 
2/51 


J/41 


P22^4 


P21  d& 
ybiP3ld6 


(5.10) 


As  in  Example  4.2  the  optimizing  D  and  the  absolute  values  of  the  elements  of 
X\  and  Y[  are 


D 


di  =  1 

'  0.427  ' 

'  0.427 

d2  =  0.82 

0.527 

0.527 

d3  =  0.41 

0.187 

0.187 

d4  =  0.63 

,     \Xi\  = 

0.269 

,     \Yi\  = 

0.269 

d5  =  0.87 

0.373 

0.373 

d6  =  0.69 

0.442 

0.279 

d7  =  0.69 

0.313 

0.465 
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with  fi(Mb)  =  6.5.  Applying  the  MPDA  property  for  block  structures  along  with  the 
relationships  of  Equation  5.10  gives  a  simple  expression  for  d6  (and  d7)  in  terms  of 
the  other  di>s  and  the  elements  of  P  of  the  form 


4  = 


d\P22 
.diP2i  . 

^ 

'd-lp33 
AlP31. 

' 

d\Pi2 
.d%Pn  . 

2 
+ 

dlPl3 

.d$Pn. 

2 

=  0.2267. 


(5.11) 


Knowing  that  d6  is  dependent  on  the  elements  of  P  and  the  other  di>s  allows  the  num- 
ber of  optimization  variables  required  for  this  problem  to  be  reduced  from  5  to  only 
4.  Variable  reduction  for  other  block  structured  uncertainties  may  be  determined  in 
an  analogous  manner  and  will  depend  on  the  number  and  size  of  individual  blocks. 
Although  Equation  5.11  is  more  complicated  than  the  relationships  for  scalar  uncer- 
tainties, the  elimination  of  unnecessary  free  variables  should  produce  a  net  savings  in 
floating  point  operations:  particularly  for  systems  of  high  order. 


5.4     Complete  Solution  to  the  Block  2x2  Problem 

In  the  original  paper  introducing  the  structured  singular  value  concept,  Doyle 
proved  that  for  k  <  3,  where  k  denotes  the  number  of  uncertainty  blocks  in  A,  the 
inequality  of 


fi(Mb)  =  sup p(MbUb)  <  m{a(DbMbDh-r) 

Ub  Db 


(5.12) 


is  guaranteed  to  attain  equality  at  the  "inf"  regardless  of  whether  or  not  a  is  sim- 
ple [15].  For  the  special  2x2  scalar  case  where  A  has  the  form 


A  = 


Si    82 
63    64 
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Daniel,  Kouvaritakis  and  Latchman  have  shown  that  equality  of  Equation  5.12  may 

still  be  guaranteed  even  though  k  =  4  [31].  The  proof  for  this  result  depends  on 
formulating  the  problem  in  a  nonsimilarity  scaling  context  and  unfortunately  does 
not  extend  directly  to  handle  block  structured  systems.  However,  by  applying  the 
reduced  scaling  structure  of  the  previous  sections  to  the  original  proof  by  Doyle,  it 
is  shown  that  the  inequality  of  Equation  5.12  always  becomes  equality  for  any  block 
structured  uncertainty  of  the  form 


A 


Ai    A2 
A3    A4 

Doyle's  original  proof  guaranteeing  equality  of  Equation  3.22  for  k  <  3  requires 
only  that  the  number  of  independent  optimization  variables  be  less  than  or  equal 
to  two  [15].  Since  one  of  the  elements  of  scaling  matrix  D  may  always  be  fixed  at 
some  arbitrary  value  without  loss  of  generality,  an  uncertainty  structure  with  three 
or  less  blocks  requires  no  more  than  two  independent  variables  for  the  optimization 
of  mioa^DMaD-1).  By  proving  that  a  2  x  2  block  structured  uncertainty  actually 
requires  no  more  than  two  optimization  variables,  this  result  may  be  extended  to 
guarantee  that  (i(Mb)  =  infi)b ^(Z^Mf,!)^1)  for  general  2x2  block  structured  systems. 

The  following  discussion  briefly  outlines  the  work  of  Doyle  leading  up  to  his  lemma 
regarding  the  equality  of  Equation  5.12  for  two  or  less  optimization  variables.  This 
outline  is  presented  only  as  a  summary  and  the  reader  is  referred  to  [15]  for  a  complete 
treatment  of  this  material. 

Before  reviewing  Doyle's  work,  a  number  of  terms  must  first  be  denned  start- 
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ing  with  the  notion  of  a  directional  derivative  for  a{DbMbDbl)  with  respect  to  the 
elements  of  D.  Denote  the  multiplicity  of  a  as  q  and  decompose  {DbMbDbl)  as 

{DbMbDb~l)  =  a(DbMbDb-')XaYaH  +  XhHhYhH 
where  Xa,  Ya  £  CnX*  for  Mb  €  CnXn.  Next,  define  some  Hj  such  that 

Hi  =  Hf  =  Re(Xf(JD6M6i)6-1)ya). 
A  directional  derivative,  V2,  may  then  be  written  as 

V2  =  {x  €  $T,  xj  =  vhHjv  |  v  e  C\  vHv  =  1} 

where  r  represents  the  number  of  elements  in  the  directional  derivative  and  similarly 
the  number  of  optimization  variables.  Formulated  in  this  manner,  V2  simply  reduces 
to  an  ordinary  gradient  for  q  =  1.  Also,  for  V2  =  0  regardless  of  the  multiplicity  of 
a,  a  true  stationary  point  is  achieved.  From  the  MPDA  theory  this  stationary  point 
assures  that  fi(Mb)  =  miD  a(DbMbDbl). 

The  convex  hull  of  V2,  denoted  coV2,  also  provides  gradient-type  information  and 
for  q  =  1,  coV2  likewise  reduces  to  an  ordinary  gradient.  However,  for  q  >  1,  a  value 
of  coV2  =  0  indicates  only  that  the  "inf"  of  infDfc  cf(A,M>iV)  has  been  achieved 
with  no  guarantee  of  a  corresponding  stationary  point.  Such  a  situation  may  actually 
be  a  "cusp"  where  fi(Mb)  <  miDb{DbMhDbl). 

Now  let  H3  =  Hf  E  Ox\  j  =  1, . . . ,  r  and  define  some  /  :  C  ->  &r 

fj(x)^xHHjx     fovxeCg. 
Also,  for  some  positive  integer  m,  let 

pm  =  {xecm\  xHx  =  i}cr 


as 
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and 


sm  =  {x  e  ar+1 1  xTx  =  1}  c  3r+1. 

Using  these  definitions,  the  following  lemma  shows  that  in  certain  cases,  V2  =coV2. 
In  such  cases,  the  condition  coV2  =  0  implies  that  V2  =  0  and  therefore  signifies  that 
a  stationary  point  has  been  found. 

Lemma  5.1:  For  /  :  C  ->  W  defined  as  above,  if  r  =  1  or  2,  then  /(P«)  is  convex. 

Proof:  See  [15]. 

Lemma  5.1  implies  that  regardless  of  the  multiplicity  of  a,  if  no  more  that  2 
optimization  variables  are  required  to  solve  Equation  5.12,  then  V2  =coV2  and  the 
minimum  of  the  infimization  must  be  a  stationary  point  and  hence  equals  (i(Mb). 
The  following  theorem  reveals  that  no  more  than  2  free  optimization  variables  are 
required  to  solve  Equation  5.12  for  systems  with  2x2  block  structured  uncertainties 
allowing  the  results  of  Lemma  5.1  to  be  invoked. 

Theorem  5.3:  For  the  case  of  2  x  2  block  structured  uncertainties,  no  more  than 
two  optimization  variables  are  required  to  solve  miDba(DbMbDb~1)  regardless  of  the 
individual  block  sizes. 

Proof:  Define  M  €  CnXn  and  2  x  2  uncertainty  matrix  A  with  corresponding  P 
matrix  as 


A  = 


Ai    A2 
A3    A4 


P  = 


P\  '  Im      Pz  "  In 


P3  "  I m      Pa  -  Im 


with  P  E  $nxn,  Pi  e  ^mxm  and  m  =  n/2.  Then  similarity  scaling  matrix  D  has  the 
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structure 


D  =  diag[du  . . . ,  <*i,  rf2j. . . ,  dz.  d3,...,d3,  d4, . . . ,  dA. 

y v '     * v '     > v '     "" v ' 

m  m  m  m 

Expanding  the  system  to  the  form  of  Equations  5.2  and  5.4  and  comparing  terms 
requires  the  elements  of  X\  and  Y\  to  be 


Xx  = 


X  771,1 

x(m+l),l 


x2m,l 


*l.lf 


r     1-23. 
x{m+l),\± 


X3m,l4 


Y,= 


2/i,i 


ym,i 


y(2m+l),l 


2/(3m),l 

u{Pi)dz 
y(2m+l),l=fp-^-i 


.  i/(3m)-1ff(p3)d4 


(5.13) 


Applying  the  MPDA  requirement  for  block  structured  uncertainties  to  Equation  5.13 
gives  a  relationship  between  the  elements  of  D  as 

d2d3 


d4  = 


a{PiP{P*) 


di   \a(P2p(P3)' 

Fixing  di  =  1  leaves  only  d2  and  d3  as  independent  variables  completing  the  proof.  ■ 
Example  5.4:  As  a  simple  example,  consider  the  following  4x4  problem  with  four 
2x2  block  structured  uncertainties. 


M  = 


with 


0.8  +  j0.3  l.O  +  jO.9  0.1+j0.5  0.4  +  ;0.5 

l.O  +  jO.4  0.7  +  j0.9  0.6  +  ;0.5  0.8  +  j0.3 

0.4  +  j0.2  0.8  +  jO.l  0.9  +  j0.3  0.5  +  jO.l 

L  0.2  + j0.5  0.7  +  j0.9  0.3  +  jl.O  0.2+j0.9 


A  = 


'  1    0 

2    0 

Ax    A2 
A3    A4 

-*  P  = 

0    1 
3    0 

0    2 
4    0 

0    3 

0    4 
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The  singular  vector  relationships  of  Equation  5.13  become 


*i  = 


xn 


321 


^31 


X41 


Xud1 


^21" 


'31* 


I   x41 


b 


>'i 


For  this  particular  example,  the  optimal  D  is 


D  = 


dl  =  1 
rfi  =  1 
c?2  =  1.38 
c?2  =  1.38 
rf3  =  1-65 
d3  =  1.65 


2/11 


2/21 


yila(F2)di 

2/21  f(TO 


2/51 


2/61 


y61nr3)d4 


d<=^=d-tJm 


d4  =  1.86  =  &&■ 


iiiiil 

(2)(3) 


with  /x(M6)  =  16.43. 


5-5     Convergence  Properties  of  the  Reduced  Scaling  Structure 

The  previous  sections  show  that  when  p,  is  found,  the  corresponding  optimal 
scaling  matrix  D  exhibits  a  degree  of  redundancy  which  allows  a  reduction  in  the 
number  of  independent  optimization  variables.  It  is,  however,  important  to  ensure 
that  when  the  relationships  between  the  dp,  and  Pij,s  are  incorporated  in  a  reduced 
scaling  structure  from  the  outset,  the  convexity  properties  of  the  original  problem 
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are  not  lost.  The  remainder  of  this  section  will  show  that,  in  fact,  these  important 
convexity  properties  do  still  apply.  Furthermore,  MPDA  conditions  used  to  establish 
equality  of  Equation  3.22  for  the  full  scaling  structure  also  pertain  to  the  reduced 
scaling  structure  guaranteeing  its  convergence  to  //. 

Convexity  of  a2(DMaD~l)  in  terms  of  eD,  D  e  V  was  first  proven  by  Safonov  and 
Doyle  [28].  Subsequent  authors  have  offered  alternative  proofs  in  somewhat  simpler 
forms  [35,  29].  The  following  lemma  concerning  the  convexity  of  some  function  / 
was  proven  by  Tsing  [36]. 

Lemma  5.2:  Suppose  Vx,  gx:  3J  i-»  3?  twice  differentiable  such  that  f(x)  =  gx(x), 
f(t)  >  9x(t)  and  -^gx(t)\t=x  >  0.  Then  /  is  convex. 

Proof:  See    [36]. 

Lemma  5.2  may  then  be  applied  to  the  scaled  singular  value  problem  to  prove 
convexity  of  a{eDMae~D).  For  notational  brevity,  denote  Mx  =  eDxMae~Dx  with 
singular  vectors  u  and  v  such  that  f(x)  =  a(Mx)  =  uHMxv.  Next,  define  gx(t)  = 
${uHeDtMae-Dtv}.  Since  f(t)  >  gx(t)  and 

d2 
-jp9S)\t=x    =    ${uH(D2Mx-2DMxD  +  MxD2)v} 

=    f{x)(uHD2u  +  vHD2v)  -  2${uHDMxDv} 


=    [uHD  vHD] 


f(x)I     -Mx 
-M*    f{x)I 


Du 
Dv 


>0. 


Therefore,  by  Lemma  5.2  /  is  convex. 

This  convexity  property  guarantees  that  any  local  minimum  of  a(DMaD-1)  with 
respect  to  D  is  always  a  global  minimum.  For  a  new  scaling  matrix,  S,  that  inherently 
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incorporates  the  reduced  scaling  structure  of  Equation  5.8,  the  convexity  lemma  still 
applies  since  S  is  a  subset  of  T>.  Therefore,  it  must  be  shown  that  the  value  of 
^(SMaS-1)  at  the  minimum  always  corresponds  to  fi  and  not  some  other  value. 

Theorem  5.4:  For  diagonal  scaling  matrices  D,  S  €  3ft™  Xn  where  S  conforms  to 
the  scaling  structure  of  Equation  5.8,  the  following  equality  must  hold 

inf  a{DMaD~^)  =  mfaiSMaS'1) 

whenever  a  stationary  point  occurs  at  the  "inf"  of  either  side  of  the  equality. 
Proof:  This  proof  requires  examination  of  the  stationary  points  at  which 
as"5 —  ~  0-  Any  maxima  of  a(SMaS~l)  with  respect  to  5  occurs  at  infinity 
because  a  increases  without  bound  with  the  free  elements  of  S.  This  eliminates 
the  possibility  of  any  stationary  point  corresponding  to  a  maxima.  Inflection  points 
are  not  possible  due  to  the  convexity  properties  of  es  so  any  stationary  point  must 
coincide  with  a  minima.  Again,  the  convexity  property  requires  any  local  minimum 
to  be  a  global  minimum  so  if  a  stationary  point  can  be  found,  it  must  correspond  to 
the  unique  minimum. 

Because  the  maximum  singular  value  is  always  a  positive  function  of  some  complex 
matrix  A  €  Cnxm,  the  stationary  points  oia2(A)  correspond  to  those  oia(A).  This  al- 
lows an  alternative  form  of  the  singular  values  to  be  written  in  terms  of  the  eigenvalues 
giving  ~cf2(A)  =  \i{AH A)  where  Ax  represents  the  eigenvalue  of  maximum  magnitude. 
Denoting  the  eigenvector  corresponding  to  Ai  as  Wi,  the  eigenvalue/eigenvector  rela- 
tionship may  be  written  as 

AHAW1  =  AiWi  (5.14) 
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and  the  derivative  with  respect  to  some  variable  x  as 

— — Wi  +  AHA— —  =  X1— —  +  W1  —  . 
ox  ox  ox  ox 

Multiplying  on  the  left  by  IVf*  gives 

ax  #x  ax        ax 

which  can  then  be  simplified  to 

giving  an  expression  for  the  derivative  of  Ai  with  respect  to  x.  The  derivative  with 
respect  to  some  element  of  scaling  matrix  S  can  then  be  found  by  replacing  the 
variables  of  Equation  5.14  such  that 

A    =    AHA=S~1M?D2MaS-1 
W1    =    Yi 
Ax    -    a2. 

The  differentiation  will  first  be  carried  out  for  a  2  x  2  case  with  respect  to  element 
522  of  S  where 


S  =  diag[slt,  s22,  s33,  S22&33 \  ]■ 

First  decompose  S  such  that 

S  =  Sq  +  S22E22  +  ■S22-533-£'44"; 


where  5*0  =  diag[su,    0,    333,    0],  k  =    /PuPii }  and  the  En  are  n  x  n  matrices  with  a 
1  at  the  z'z'th  position  and  0  elsewhere.  Likewise,  S'1  may  be  decomposed  as 

~       °  «  h    •       «       fr" 

-522  •S22'S33/>' 


Differentiation  of  A  with  respect  to  s22  then  proceeds  directly  with 
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dA 
ds 


22 


■^Mf&MaS-1 


E. 


44        l\/rH  c2 


'22 


-S~1M"S2Ml 


a    *-/    '"o 


■s22,S33fc 
^22  £*44 


M?S2MaS- 


'22  ,322 


s22533"; 


+2522S,-1Mf  feikf^"1  +  EAisl3k2MaS~l 


The  jE,-,-  terms  contain  only  one  nonzero  value  so  this  expression  may  be  rewritten  as 


8A 
ds 


•S22 


22 


=    -^E22S-1M^S2MaS~l- 


"22 


•S22-S33 
522533"; 


-1  fl/f#  C2  n/T    C-l 


E44S~1M"S2MaS 


522C-lwffc2 


sii  si2s33k 


+  2s22S-lM^E22MaS-1  +  2<S224^2S-1Mf  E^MaS-\ 


Multiplying  ^  on  the  left  by  Yf  and  the  right  by  Y\  as  in  Equation  5.15,  gives  an 
expression  for  J^-  of 


«9^  =  y„  ^4 

<9<S22  *     5^22 


yH  yH 


-2 

^22  -S22 

+2s22Y1H  \s-lM*E22MaS-1  +  s^S-xM? EuMaS~l\  Yt. 


Noting  that 


and 


—  -T-— 

ds22  ds22 


rH 


Y?S^MZSE%i  =  oX?El% 


a  final  expression  for  ■—  may  be  written  as 


3*22 


da  a 

ds22      s22 


X?(E22  +  E,4)X1  -  YXH{E22  +  E44)YX 


(5.16) 
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Carrying  out  this  process  with  respect  to  s33  gives  an  expression  similar  to  Equation 

5.16  of 

da         a 


ds33       s33 


X?  (£33  +  E44)X,  -  YfiEss  +  E44)Y1\  .  (5.17) 


Setting  Equations  5.16  and  5.17  to  zero  establishes  conditions  at  which  the  new  scaling 
technique  reaches  its  global  minimum.  The  stationary  point  therefore  corresponds  to 
the  value  of  S  where 

[^21^21  -  y21y2i  +  x41x41  -  y4ly4r}  =  0  (5.18) 

and 

[X3l2>31  -F31S/31  +^41^41  -  y4iJ/4i]  =  0  (5.19) 

simultaneously. 

The  MPDA  property  guarantees  that  when  element -by- element  alignment  of  the 
principal  singular  vectors  of  o=(5,MaS'-1)  occurs,  some  Ud  exists  such  that  a(5'Afa5'-1)  = 
p(MaUd)  =  ft(Ma).  From  Equations  5.18  and  5.19,  such  an  alignment  of  X\  and  Yx 
always  coincides  with  a  stationary  point  of  ^{SMaS-1)  because  it  requires  that 

M2  -  I2/21I2  =  |x3i|2  -  |y3i[2  =  |^4i|2  -  |2/4i|2  =  0. 

Therefore  one,  and  hence  the  only,  minima  of  a(SMaS~l)  with  respect  to  S  must 
correspond  to  fi(Ma)  as  long  as  some  stationary  point  does  exist. 

Extension  to  higher  order  systems  is  simply  a  matter  of  determining  which  ele- 
ments of  S  depend  on  a  particular  s{.  For  each  occurrence,  an  additional 

-%i  EjjXi  —  Yx  EjjYj 
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expression  is  included  where  the  jj  terms  denote  the  elements  of  S  dependent  on  a 
particular  ${  term.  Again,  general  n  x  n  systems  similarly  invoke  MPDA  conditions 
guaranteeing  that  the  stationary  point  corresponds  to  n(Ma).  ■ 


CHAPTER  6 
DIRECT  RELATIONSHIP  BETWEEN  R,  L,  AND  D 

6.1     Problem  Formulation 

The  results  of  Chapter  5  show  the  number  of  optimization  variables  required  for 
similarity  scaling  to  be  no  more  than  that  required  for  nonsimilarity  scaling.  It  would 
therefore  seem  consistent  that  some  direct  relationship  exists  between  scaling  matrices 
D  for  mtDtvaiDMaD-1)  and  R  and  L  for  infi?)La(ir1ML-1)ff(£P.R)  since  both 
methods  accurately  compute  fi.  The  following  Theorem  presents  such  a  relationship. 
Theorem  6.1:  For  element-by-element  structured  uncertainties,  given  some  scaling 
matrix  D  conforming  to  the  structure  of  Equation  5.8,  a  corresponding  R  and  L  of 
the  form 


and 


dt  +  dn+i  +  d2n+i  •  •  •  an(n_1)+t- 


Pll       I      P212  PJn 

$  z d?     < 

P?,l  I  Ph  Pin 


n 


(6.2) 


tP  '       H2  "  ri'i 

n(i-l)+l  an(i-l)+2  an(i-l)  +  n 

may  be  found  without  further  optimization  such  that 


*i{R-1ML-1)a(LPR)  =  ar(DMaD-1)    i  =  l,...,n.  (6.3) 

Proof:  Rather  than  dealing  explicitly  with  singular  values,  it  is  more  convenient 
in  this  case  to  use  the  equivalent  representation 

Vi{A)  =  ^\i{A"A)    AeCnXn,    i  =  l,...,n 
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and  work  with  the  characteristic  equations  of  A  =  {R~l  M  L'1)11  {R~l  M  L~l)  and 
$  =  {DMaD-1)H(DMaD~l).  The  remainder  of  this  proof  will  consider  the  2  x  2  case 
for  notational  simplicity  but  the  results  will  be  shown  to  hold  for  the  general  n  x  n 
case  as  well.  Define  R  =  diag[l,  r2],  L  =  diag[l,  l2],  D  =  diag[dt,  d2,  d3,  d4],  and 


M  = 


11 

M12 

21 

M22 

N   Ifcl 


<  P. 


Matrix  Ma  is  formed  by  expanding  M  and  P  such  that  the  elements  of  P  are  diago- 
nalized.  Ma  therefore  retains  the  rank  of  M  so  the  characteristic  equations  of  A  and 
/3  may  be  written  respectively  as 

A   -*      X2  +  bX  +  c 
P    ->      X2  +  BX  +  C. 


In  order  for  Equation  6.3  to  hold  for  all  i  simultaneously,  it  is  necessary  that 
bcf2(LPR)  =  B  and  ca4(LPR)  =  C.    Expanding  A  while  retaining  the  individual 
elements  of  i?,  L,  and  M  gives 


A  = 


Afn  Mn  +  M2iM2i     MuMn.  j.  M21M21 


2 '2 


M12MV    J_  M22M21        M12M,2     1    M22M22 


so  that 


MuMn  +  MnM21  4.  Mi2Mi2      M22M22 


/I 


„2/2 
'  2*2 


and 


c  = 


_  M11M11M22M22  +  M21M21M12M12  -  MUM12M22M21  -  MXXMX2M2XM; 


22 


r.2/2 

2*2 
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By  similarly  expanding  /5,  B  and  C  are  found  to  be 


B 


Mil Afn  [di  +  d§ 
Mi2M12  [d2  +  d23 


,  P12" 
4. 

p1i 
«8 

d| 

-M21M21  4  +  d 


M22M22  d2  +  d 


rii 

+ 

pL" 

<*? 

c3. 

P21 

+ 

P22' 

«§ 

dt 

C    =    K-   MnMnM22M22  +  M21M21M12M12 


K  ■   M11M12M22M21  +  MiiMi2M2iM22 


A"    = 


d?  +  ^    •   d2  +  d; 


P11P21    ,  P11P22  _>_  P12P21    ,  P12P22 

d\d\     did\  +  d^i  +  d\d\ 


Equating  ba2(LPR)   =  B  and  ca4(LPR)   -  C,  the  following  direct  relationships 

between  R,  L,  and  D  will  hold 

2  _l  ^21  [4+4?" 

(6.4) 


r2  - 

T2  ~ 

~dl  +  d23~ 
d\  +  dl 

I2  - 

1            l2  ~~ 

P11    1   P12 

d\    T    <** 

p21       I      P22 

.  d§    +    dl  _ 

but  only  when  equality  constraint 

a\L 

PR)  = 

dl  +  dl]  ■ 

Pll     ,    P12" 

(6.5) 

applies.  At  first  glance  it  may  appear  that  the  constraint  of  Equation  6.5  severely 
limits  conditions  under  which  Equation  6.4  holds.  The  remainder  of  the  proof  will 
reveal,  however,  that  this  constraint  simply  requires  the  reduced  similarity  scaling 
structure  developed  in  Chapter  5. 

The  left  hand  side  of  Equation  6.5  may  always  be  represented  as  Ax  [(LPR)H(LPR)]. 
Expanding  matrix  (LPR)H(LPR)  gives 


(LPR)H(LPR)  = 


Pll  +  Pll12  PUP12T2  +  P2lP22%T2 

P\\Pi2T2  +  P21P22/2S         Pl2i%  +  v\2r22l\ 
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The  corresponding  characteristic  equation  of  this  matrix  becomes 


y  fcb\  +  c.  - 


Denote  the  two  roots  of  Equation  6.7  as  a2(LPR)  and  a2{LPR).  Coefficients  6  and 
c  may  be  expressed  in  terms  of  these  roots  such  that 


b    =    a2(LPR)  +  a2(LPR) 
c    =    a2(LPR)  ■  a2(LPR). 


(6.6) 
(6.7) 


Relating  these  expressions  to  the  elements  of  (LPR)  gives 


b  =  a\LPR)  +  aj(LPR)  =  p2n  +  p221l2  +  p\2r\  +  v\2r\l\ 


(6.8) 


and 


c    =    a2(LPR)  ■  a2(LPR) 

=    (Pn  +  p\A)  ■  (Pn4  +  p222r22l2) 

~  (PllPl2r2  +P2lP22llr2)  ■  {puPi2r2  +  P21P22&2) . 


(6.9) 
(6.10) 
(6.11) 


Replacing  r2  and  l2  in  Equations  6.8  and  6.11  with  the  relationships  of  Equation  6.4 
and  simplifying  gives 


b  =  fix  +  Pk 


d2  +  d2 


and 


1^2  +  <%} 


di  +  4 


+  Pn 


jg     +    g 

P21       I      P22 


+  P22 


4  +  ^ 


d\     "*"     rf| 
^21       I      P22 


Ml      .      PL 


[PllP22    -  Pl2P2l]" 


[^  +  <8]  • 


p|i  _i_  P22 

I    +    d? 


(6.12) 
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Equation  6.7  shows  coefficient  c  to  be  the  product  of  a2(LPR)  and  a2{LPR).  Com- 
bining Equations  6.5  and  6.12  then  gives  the  expression  for  cr$(LPR)  of 


\pnP22  ~  P12P21] 


d\    ^    ft 


(6.13) 


[d2  +  dj] 

Equation  6.7  also  shows  that  coefficient  b  is  the  sum  of  a2(LPR)  and  a2(LPR) 
so  that  b  may  likewise  be  written  as 


b  = 


d\  +  4 


Pu    ,  P12 
d\  +  <*g 


+ 


JP11P22  -  P12P21Y 


"22 n  ' 

24 


This  gives  a  new  expression  for  6  that  must  be  satisfied  along  with  Equation  6.8 

that 


so 


b    =    Pu+Pu 

=    d\  +  d\ 


d\  +  d\ 


dl  +  dl 


P11    ,   P12 
[d\  +  d2\ 


£11   1   P12 
j%  +  jg 

P21       I      P22 

W  +  if 


P21  1  P22 


K2  +  <S] 

Multiplying  out  these  two  expressions  and  eliminating  common  terms  reveals  that 
they  are  equal  whenever 


d\d\p\2p221  +  d\d\p\xp\2  -  2d2d22dld24puPl2p21p22  =  0. 


(6.14) 


Equation  6.14  can  be  rewritten  as 


]2 


did4Pi2P2i  -  d2d23pnp22 


so  the  constraint  of  Equation  6.5  will  hold  for  D  and  P  such  that 


d4  =  d2d3. 


/P11P22 
P12P21 
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This,  however,  is  the  same  dependent  relationship  between  the  pijis  and  diis  of  Equa- 
tion 5.8  so  the  direct  relationship  between  R  and  L  of  nonsimilarity  scaling  and  D  of 
similarity  scaling  will  always  exist  whenever  the  new  scaling  structure  is  invoked. 

For  the  n  x  n  case,  a  general  expression  for  R  and  I  as  a  function  of  D  may 
be  derived  in  an  analogous  manner  by  comparing  the  characteristic  equations  of  the 
corresponding  matrix  systems.  For  M  €  C3x3,  the  characteristic  equations  would 
contain  three  coefficients  which  must  then  be  matched  up  as  before.  While  intuitively 
simple,  the  actual  equations  quickly  become  too  unwieldy  to  present  here.  However, 
a  general  expression  similar  to  that  of  Equation  6.4  may  be  written  as 


,2  _   "l  +  "w+1  +  d2n+l  •  •  •  dn(n-i)+i        .  _ 

! '  d2  -I-  d2       4-  A2         .  .  .  A2  I  —  L  .  .  . 

"»   T  Un+j  "I"  «2n+i  "n(n-l)+t 


n 


and 


/?  = 


Pll       1      Pl2 

<q  +  afij 

Pin 

P?,l              1             P?,2 
d2                     '     rf2 

an(i-l)+l              an(i-l)  +  2 

„2 

n(t  — l)  +  n 

l...n 


completing  the  proof.  ■ 

The  following  example  illustrates  the  application  of  these  expressions. 

Example  6.1:  Using  M  and  D  from  Example  4.1,  optimal  values  for  R  and  L  are 
found  to  be 


r2 


d\  +  d\  +  d? 

Vl  +  ^l  +  d| 


=  1.151, 


rz 


dj  +  dj  +  <3P7 
^dl  +  dl  +  dj 


=  1.174 


and 


\ 


Pll       I      Pl2      I      Pl3 
J2        I       J2     T      j2 


3? 


<S 


3 


P21       I      P22      1      P23 
2      T"     J2     T     J2 


=  0.9443,    I3  = 


3 


3 


dl 


\ 


Pn 
d\ 


I      Pl2      I      Pl3 


^31       I      P32      I      P33 


-  1.203 
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with 


a^R-'ML-^aiLPR)    =    ^{DMaD~l)  =  n(Ma)  =  8.25 
ai{R-1ML-x)a{LPR)    =    a2{DMaD-1)  =  5.72 
a^R-^ML-'piLPR)    =    az{DMaD-1)  =  1.76. 

It  should  be  noted  that  this  transformation  not  only  provides  L  and  R  given  D, 
but  also  allows  the  computation  of  D  given  R  and  L.  However,  this  reverse  trans- 
formation results  in  a  transcendental  function  that  must  be  solved  iteratively.  The 
iterative  solution  is  still  quite  fast  since  it  involves  only  elementary  math  operations: 
no  singular  value  or  eigenvalue  decompositions. 

6.2     Independence  from  M 

An  interesting  observation  concerning  Equations  6.1  and  6.2  is  that  they  are  not 
functions  of  M.  This  independence  from  M  requires  that  any  P,  R,  L,  D  combination 
satisfying  Equation  6.3  for  a  given  M,  must  also  satisfy  Equation  6.3  for  any  other 
M.  As  a  simple  example  consider  the  following  system  with  two  randomly  generated 
M  matrices. 


Afi  = 


0.87  -  £2.01      4.70  -  £0.08 


4.16  — £0.59    —0.92  +  £0.13 


r                -i 

M2  = 

4    3 

2    1 

p 


1.59    1.64 


D  =  diag[l,  1.415,  1.401,  1.057] 


2.81    0.83 
where  D  satisfies  the  new  scaling  structure.  Values  for  R  and  L  are  found  to  be 

R  =  diag[l,  0.974]    L  =  diag[\,  0.912] 
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so  that  Equation  6.3  becomes 


<yi(BrlM1L-1p{LPR)    =    a1{DMalD~l)  =  20.28 
<x2(R-1M1L-1p(LPR)    =    a2{DMalD-1)  =  UM. 

Using  the  same  R,  L,  and  D  values  for  M2  with  no  optimization  gives 

cr1(R-1M2L-1)a{LPR)    =    ^(M^i)-1)  =  19.18 
^(R-iMiL^piLPR)    =    a2{DMa2D-l)  =  l.M. 

The  direct  relationship  between  R,  L  and  £>  leads  to  the  following  corollary  con- 
cerning systems  that  "cusp"  during  the  optimization  process. 

Corollary  6.1:  Systems  with  element-by-element  structured  uncertainties  with  non- 
^}V}e°{R~~lML-1)  ("cusps")  for  nonsimilarity  scaling  must  also  produce  cusps  for 
similarity  scaling  techniques  and  vice-versa. 

Proof:  The  proof  follows  directly  from  the  one-to-one  correspondence  between  the 
ffiu  oicri{R-lML-l)a(LPR)  and  ai{DMaD-x)  required  by  Equation  6.3.  ■ 

Corollary  6.1  serves  to  eliminate  the  appealing  notion  that  a  one  of  these  two 
scaling  techniques  may  converge  to  ft  even  if  the  other  method  fails.  While  some 
alternative  scaling  technique  may  eventually  overcome  the  difficulties  associated  with 
cusping  systems  it  is  the  case  for  nonsimilarity  or  similarity  scaling. 

Besides  ruling  out  the  use  of  an  alternative  scaling  method  for  cusping  systems, 
the  direct  relationship  between  R,  L,  and  D  serves  to  improve  the  vector  optimization 
method  of  Fan  and  Tits.  The  next  chapter  includes  a  discussion  of  this  method  and 
the  improvements  offered  by  the  R,  L,  D  relationship. 


CHAPTER  7 
REDUCTION  OF  OPTIMIZATION  VARIABLES  REQUIRED  FOR  THE 

METHOD  OF  FAN  AND  TITS 

7.1     Constrained  Vector  Optimization 

While  the  scaling  structure  developed  in  Chapter  5  allows  a  reduction  in  the  num- 
ber of  optimization  variables  required  to  compute  /i,  the  process  still  involves  a  large 
number  of  singular  value  decompositions.  To  avoid  these  computationally  intensive 
decompositions,  Fan  and  Tits  introduced  a  constrained  vector  optimization  proce- 
dure that  calculates  \i  using  a  series  of  vector  norms  [37].  The  problem  formulation 
is  similar  to  the  eigenvalue/eigenvector  relationship  for  matrix  A  £  CnXn 

AWi  =  XtWz 

where  W{  represents  an  eigenvector  and  A;  the  corresponding  eigenvalue.  Since  fj,(Ma) 
is  actually  the  solution  to  supw  p(MaU),  there  exists  some  supremizing  U  such  that 

||MaC7Wi||2  -  ||A1(Ma/7)W1||2  =  P{MaU)\\Wih  =  KMa). 

For  ||Wi|[2  =  1,  this  reduces  to  UilfB#Wi||a  =  p{MaU)  =  fi(Ma).    Defining  vector 
x  —  UW\  then  gives 

\\Max\\2  =  p:(Ma).  (7.1) 

Computing  the  2-norm  of  a  complex  16  x  16  vector  requires  about  4n  or  64  flops  versus 
more  than  50,  000  flops  for  the  singular  value  decomposition  of  a  complex  16  x  16 
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matrix.  Therefore,  optimizing  the  left-hand  side  of  Equation  7.1  offers  potentially 
significant  numerical  advantages  in  the  calculation  of  fi. 

The  optimization  must  first  be  constrained  to  include  only  "2"  vectors  of  the  form 
UW(.  This  constraint  is  imposed  by  the  Fan  and  Tits  formulation  for  (s(Ma): 

fi(Ma)  =  max  {||Max||a  |  ||P;x||2||Maz||2  =  \\ViMax\\2,  i  =  l,...,m}  (7.2) 

x£C"2 

where  the  Vi*s  are  projection  matrices  with  block  size  m  and  x  G  C"2.  Unfortunately, 
this  function  is  non convex  so  it  offers  only  candidate  solutions.  Determining  whether 
these  candidates  equal  fi(Ma)  involves  computing  a  corresponding  scaling  matrix  D 
from  x  and  ensuring  that  ||Max|j2  =  a(DMaD-1).  Also,  in  its  original  formulation, 
this  expression  requires  n2  -  1  optimization  variables.  The  following  Theorem  shows 
that  the  scaling  structure  introduced  in  Chapter  5  applies  to  Equation  7.2  so  the  Fan 
and  Tits  method  actually  requires  no  more  than  2(n  —  1)  optimization  variables. 

Theorem  7.1:  Given  M  €  Cnxn  and  A  €  Cnxn,  where  A  is  composed  of  1  x  1 
blocks,  the  solution  for  /i(Ma)  in  the  constrained  optimization  problem 

n{Ma)  =  max  {||Maz||2  |  ||P,-x||2||Maa:||2  =  \\ViMax\\2,  i  =  1, . . .  ,m} 

with  simple  a  can  be  found  with  no  more  than  2(n  —  1)  free  variables. 

Proof:  As  in  Section  5.5,  denote  the  reduced  scaling  structure  of  Equation  5.8  as 
S.  Let  S  be  the  solution  to  misa(SMaS'1)  and  Yx  the  corresponding  right  singular 
vector.  Then  the  optimal  solution,  x,  to  Equation  7.2  may  be  written  as    [38] 

5-^1 


x  — 


\\s-wA 


12 
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The  dependent  relationships  between  the  elements  of  Yx  may  be  determined  using 

the  method  of  page  46.  Since  Yx  is  a  singular  vector  of  SM.S'1,  any  scalar  multiple 

of  Yx  will  also  be  a  singular  vector  of  SMaS~l.  Therefore  it  can  be  scaled  such  that 

2/n  =  1  without  loss  of  generality.  Likewise,  ia  may  also  be  scaled  to  1.  Carrying  out 

the  multiplication  of  S  and  Fi,  results  in  a  general  expression  of  the  form 


ZT1^  = 


1 

<—  fixed 

2/2,1 

<—  free 

2/n,l 

'. 

Vn+1,1 

<—  free 

2/2,1^+1,1 

<—  fixed 

2/n,iyn+l,l 

<—  fixed 

2/271  +  1,1 

«—  /ree 

2/2,12/271+1,1 

<—  fixed 

(7.3) 


J/n,l2/(n-l)n+l,l    J      <—  fixed. 

Subtracting  the  (n  -  l)2  fixed  variables  of  Equation  7.3  from  the  n2  -  1  free 
variables  previously  required  leaves  the  desired  result  of  n2  -  1  -  (n  -  l)2  =  2(n  -  1) 
free  optimization  variables.  ■ 

Example  7.1  illustrates  the  application  of  Equation  7.3  to  the  system  of  Example 


4.1. 
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Example  7.1:  Using  the  system  of  Example  4.1,  the  optimal  S^Yi  may  be  written 


as 


S'^Yi  = 


Scaling  j/n  and  Si  to  1  gives 


S"1*!  = 


with  corresponding 


y^m 


y^m 


«4 


2/41-2 


Pl2 


J/41 12" 


*2PHS4 
P21 


S|pilS4 


3/71 
S7 


Vn-a 


P21 


SJPUS7 


1     ■ 

r         1 

2/21 

1.236 

2/31 

1.094 

2/41 

0.855  +  0.519; 

2/212/41 

= 

1.056  +  0.641; 

2/312/41 

0.935  +  0.567; 

2/71 

0.355  +  0.935; 

2/212/71 

0.438  +  1.155; 

2/312/71    . 

.  0.388  +  1.023; 

||M0z||2  =  8.25=Zf(Z>Af0Zr1). 


7.2     Application  of  the  Relationship  Between  R,  L  and  S 


As  mentioned  earlier,  the  solution  to  Equation  7.2  is  nonconvex  with  many  local 
maxima  possible.    The  only  means  of  ensuring  that  a  candidate  solution  vector,  x, 
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actually  corresponds  to  fi(Ma)  is  to  determine  optimal  scaling  matrix  D  such  that 

||M«a||3  =  aiDM.D-1)  =  ft(Ma).  (7.4) 

(The  method  for  directly  computing  D  given  x  is  presented  in  Section  8.6.)  Unfortu- 
nately, this  singular  value  decomposition  must  be  performed  on  the  expanded  n2  x  n2 
system  for  each  candidate  solution,  reducing  the  overall  efficiency  of  this  technique. 
However,  using  the  direct  transformation  from  D  to  R  and  L  presented  in  Chapter 
6,  the  verification  test  for  fi(Ma)  may  be  reformulated  as 

\\Max\\2  =  a(R-1ML-1)a{LPR)  (7.5) 

where  the  left-hand  side  consists  of  only  n  x  n  matrices.  For  a  system  with  n  =  3, 
the  expanded  n2  x  n2  matrices  of  Equation  7.4  require  about  10,  000  flops  for  a  single 
singular  value  decomposition.  Equation  7.5,  on  the  other  hand,  requires  about  1,000 
flops  including  both  singular  value  decompositions  and  the  required  transformation 
from  D  to  L  and  R.  Even  more  substantial  savings  occur  as  n  increases. 

The  Fan  and  Tits  method  outlined  in  this  chapter  first  appeared  in  1985  [37]. 
Their  constrained  optimization  procedure  was  implemented  using  specialized  soft- 
ware from  the  proprietary  Harwell  Subroutine  Library  [39].  Since  that  time,  the 
growing  acceptance  of  the  structured  singular  value  concept  has  led  to  a  number  of 
alternative  formulations  of  /j,.  These  formulations  are  aimed  at  improving  its  com- 
putational properties  to  allow  for  interactive  design  without  the  restrictions  imposed 
by  the  use  of  proprietary  software.  The  recent  announcement  by  the  Math  Works  of 
the  "^-Analysis  and  Synthesis  Toolbox"  for  MATLAB™  indicates  the  demand  for 
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user-friendly  software  to  calculate  p  [40].  This  new  toolbox  uses  a  variation  of  the 
power  method  proposed  by  Packard,  Fan,  and  Doyle  to  compute  a  lower  bound  for 
fi(Ma)  [41].  This  optimization  searches  for  stationary  points  of  the  singular  values  of 
(DMaD-1)  with  respect  to  the  elements  of  D.  From  the  Principal  Direction  Align- 
ment theory  of  Daniel,  Kouvaritakis,  and  Latchman,  such  stationary  points  have  been 
shown  to  coincide  with  stationary  points  of  P{MaUd)  with  respect  to  the  elements  of 
Ud  [31].  However,  like  all  known  lower  bound  expressions,  the  procedure  is  noncon- 
vex  so  candidate  solutions  for  fi(Ma)  must  be  compared  to  the  singular  value  upper 
bound  to  see  if  they  actually  achieve  the  true  supremum.  Again,  rather  than  per- 
forming a  singular  value  decomposition  on  the  n2  x  n2  system  DMaD~\  the  test  may 
be  formulated  as 

KM*)  <v(R-1ML~1)v(LPR)  (7.6) 

involving  only  n  x  n  matrices  where  jl{Ma)  represents  a  candidate  for  fi(Ma),  If  the 
inequality  of  Equation  7.6  holds  with  equality,  then  fi(Ma)  =  ft(Ma). 

If  no  candidate  for  \x  is  found  for  which  the  lower  and  upper  bounds  become  equal, 
then  the  system  may  have  a  repeated  maximum  singular  value  at  the  "inf "  of  Equation 
5.1  for  scalar  uncertainties  or  Equation  5.9  for  block  structured  uncertainties.  Lack 
of  stationarity  at  this  upper  bound  implies  that  no  p,  exists  for  which  the  upper 
and  lower  bounds  become  equal.  Chapter  8  contains  a  detailed  description  of  the 
difficulties  involved  with  computing  fx  when  the  maximum  singular  value  repeats. 


CHAPTER  8 

CALCULATION  OF  THE  STRUCTURED  SINGULAR  VALUE  WITH 

REPEATED  MAXIMUM  SINGULAR  VALUES 

8.1     Cusping  Singular  Values 

The  work  presented  in  the  previous  chapters  generally  requires  that  the  maximum 
singular  value  not  repeat.  This  requirement  stems  from  the  fact  that  the  inequality 
in  the  right-hand  side  of 

p(Ma)  =  sup p(MaUd)  <  mfaiDMaD'1) 

Ud  D 

is  not  guaranteed  to  achieve  equality  at  the  "inf"  if  a  -  <ru  i  ^  1  because  this 
suggests  that  a  stationary  point  with  respect  to  D  has  not  been  found.  Following  the 
infimization  oia{DMaD~1)  with  respect  to  D,  three  possible  conditions  can  occur. 
These  three  conditions  are  illustrated  in  Figures  8.1  through  8.3  with  corresponding 
Ma  matrices  in  Tables  8.1  through  8.3  respectively.  Figure  8.1  represents  a  system 
in  which  MD  a{DMaD~l)  is  achieved  with  distinct  a(DMaD^).  From  the  MPDA 
property  reviewed  in  Chapter  4,  this  guarantees  the  simultaneous  determination  of 
supUdp(MaUd)  and  hence  fi(Ma). 

Figure  8.2  represents  a  second  system  where  infD ^{DMaD'1)  occurs  with  a  = 
a2  =  sup^  p(MaUd).  The  point  at  which  the  equalities  hold  will  be  referred  to  as  a 
"kiss"  and  its  calculation  requires  special  consideration. 

Figure  8.3  shows  a  system  in  which  infD  a(DMaD^)  occurs  with 
a  =  <r2  ^  supUd  p(MaUd)  =  p{Ma).  This  case  will  be  referred  to  as  a  "cusp"  and 
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Figure  8.1:     Noncusping  System. 
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23.5 
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Figure  8.2:     System  with  "Kiss". 
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11.5 


Figure  8.3:      System  with  "Cusp". 
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-0.89  +  2.04; 

3.87  +  17.96.; 

-7.68  -  3.20; 

-5.11  -11.27? 

9.83  -  18.68; 


Table  8.1:  Ma  Matrix  for  Figure  8.1. 


-4.34  -  5.31; 

1.66  -  4.72; 

-11.10-8.20; 

-2.81  +  14.23; 

-5.39  +  6.77; 


-1.83  +  4.69.; 

6.83  +  2.94; 
6.84  -  20.10.7 
8.19  +  11.80; 

1.11  +  0.19; 


3.49  -  6.31; 

-7.38  -  14.37; 

-3.15  +  0.62j 

8.89  +  14.60; 

6.29  +  2.55; 


3.07  +  5.92; 
-8.73  -  1.29; 

2.27  +  7.21; 
2.73  -  10.69; 
-2.83  +  6.75; 


Table  8.2:  Ma  Matrix  for  Figure  8.2. 


19.18  +  0.37; 
-0.20  -  3.07; 

6.42  +  1.85; 

0.06  +  0.64; 
-2.39  -  5.34; 


6.82  -  1.75; 

18.56  +  1.29; 

-0.70  +  1.03; 

-0.53  -  2.47; 

3.21  -  0.78; 


3.13-0.95; 

-1.44  +  0.35; 

15.34-2.01; 

3.53  -  0.97; 

3.74  +  1.38; 


-4.92  +  1.11; 

3.22  +  2.37; 

-0.77  -  0.82; 

10.96-3.11; 

3.24  -  0.03; 


3.34- 

-  4.59;  " 

-1.32  +  3.15; 

-0.13  +  1.36; 

3.93- 

-  0.96; 

15.41  - 

-  0.33; 

5.18  +  0.37; 
-0.20  -  3.07; 

6.42  +  1.85; 

0.06  +  0.64; 
-2.39  -  5.34; 


Table  8.3:  Ma  Matrix  for  Figure  8.3. 


6.82  -  1.75; 

4.56  +  1.29; 
-0.70  +  1.03; 
-0.53  -  2.47; 

3.21  -  0.78; 


3.13  -  0.95; 
-1.44  +  0.35; 
1.34-2.01; 
3.53  -  0.97; 
3.74  +  1.38; 


-4.92  +  1.11; 

3.22  +  2.37; 
-0.77  -  0.82; 
-3.03-3.11; 

3.24  -  0.03; 


3.34  -  4.59; 
-1.32  +  3.15; 
-0.13  +  1.36; 

3.93  -  0.96; 

1.41-0.33; 


calculating  n{Ma)  for  such  a  system  represents  a  quite  challenging  problem.  The 
remainder  of  this  chapter  will  address  several  approaches  that  greatly  improve  upon 
previous  techniques  for  determining  p{Ma)  when  a{DMaD~l)  repeats. 

As  previously  discussed,  noncusping  systems  such  as  that  of  Figure  8.1  are  guaran- 
teed to  achieve  p  at  info  a(DMaD^).  Using  the  following  theorem  due  to  Latchman 
[27],  it  is  possible  to  compute  the  supremizing  Ud  of  suP[/d  p(MaUd)  given  the  infimiz- 
ing  D  of  miDa(DMaD-1). 

Theorem  8.1:  Given  the  optimal  diagonal  scaling  matrix  D  corresponding  to  a 
stationary  point,  there  exists  a  diagonal  unitary  matrix,  Ud  =  dingle3*1 , . . . ,  e^],  0  < 
9i  <  2tt,  where  0,-  =  [arg(yil)  -arg(x,i)],  i  =  1, . . . ,  n,  such  that  the  principal  singular 
vectors  Xx  and  Yx  of  DMab-1Ud  are  aligned. 

Proof:  Recall  that  p{MaUd)  <  a{DMaD-lUd).  For  any  system,  singular  vectors 
Xi  and  Yx  obey  the  following  relationships: 

{DMab-l)Yx    =    aXx 
X?0Mutr*)    =    aY".  (8.1) 

Introducing  a  diagonal  unitary  matrix,  Ud,  allows  Equations  8.1  to  be  rewritten  as 

(DMaD-'U^UfYx    =    aXx 
X?(bMab-'Ud)    =    aYfUd.  (8.2) 

From  Equations  8.1  and  8.2  it  is  apparent  that  Ud  may  alter  the  phases  of  Xt  and 
Yx  while  continuing  to  satisfy  the  MPDA  requirement  that  the  magnitudes  of  the 
corresponding  elements  of  X1  and  Yx  be  equal.    The  independent  phase  alignment 
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of  the  elements  of  Xx  and  Yt  provides  the  necessary  information  to  compute  the 
supremizing  Ud  (denoted  as  Ud).  ■ 

Using  Theorem  8.1,  the  individual  diagonal  elements  of  Ud  are  found  to  be 

for  element-by-element  structured  uncertainties.  A  similar  expression  reflecting  indi- 
vidual uncertainty  blocks  applies  to  the  general  block  structured  case.  The  actual  Ud 
determined  in  this  manner  is  unique  within  some  scalar  multiplier  e^. 

Example  8.1:  Applying  the  relationship  of  Equation  8.3  with 
D  =  diag[l,  0.460,  0.350,  0.425,  0.541],  the  corresponding  Ud  for  the  system  of  Figure 
8.1  is 

Ud  =  diag[l,  0.996 -j0.078,  -0.998 +  j0.060,  0.966  +  j0.250,  0.921  -j0.390] 

with  a{DMaD-i)  =  p(MaUd)  =  36.94  =  fx(Ma). 

The  ability  to  directly  compute  Ud  given  D  provides  a  simple  means  to  ensure  that 
aiDMaD-1)  =  n(M9)  because  it  guarantees  that  suP[/d  p(MaUd)  =  infD ^DMaD'1). 
This  becomes  increasingly  important  as  the  conditions  illustrated  in  Figure  8.2  are 
considered. 

8-2     Systems  with  infg  ajDMaP-1)  =  a2{DMaD-r)  =  p(Ma) 

When  the  singular  values  become  equal  or  nearly  equal,  the  corresponding 
singular  vectors  become  unstable  with  small  changes  in  the  ails  resulting  in  large 
variations  in  the  columns  of  X  and  Y    [42].    This  behavior  is  a  direct  result  of  the 
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expressions  defining  the  singular  vectors  of  a  matrix  A,  namely 

AAHXi    =    afXi 

AHAYi    =    a}Yz.  (8.4) 

When  the  singular  values  repeat,  the  singular  vectors  combine  to  span  a  g-dimensional 
subspace  where  q  represents  the  multiplicity  of  the  repeating  o>5.  For  q  =  2,  (the 
most  common  case),  the  first  line  of  Equation  8.4  becomes 

AAHXX     =    <j\Xx 

AAHX2    =    a\X2  (8.5) 

with  the  resulting  subspace  formed  by  the  linear  combination  of  Xx  and  X2  to  give 

AAH[(i  +  J6)X1  +  (a  +  jfi)X2]  =  <rf  [(7  +  j6)X1  +  (a  +  j/3)X2}.  (8.6) 

(Note  that  Y\  is  dependent  on  X{  through  the  relationship  AHXt  =  afc  and  thus  a 
similar  development  for  the  subspaces  spanned  by  Yx  and  Y2  would  be  redundant). 
From  Equation  8.6  it  is  apparent  that  any  linear  combination  of  Xx  and  X2  satisfies 
the  definition  of  a  singular  vector.  This  arbitrariness  poses  a  significant  difficulty  for 
algorithms  relying  on  the  vector  alignment  properties  of  MPDA  theory  because  there 
is  no  longer  any  guarantee  that  vectors  chosen  arbitrarily  from  the  ?-dimensional 
subspace  will  satisfy  the  alignment  requirements. 

For  example,  the  system  of  Figure  8.2  has  miD  u{DMaD'x)  =  24.14.  Using  the 
method  of  Theorem  8.1  for  Xu  Yx  and  X2,  Y2  gives  corresponding  values  of  p(MaUd) 
as  follows: 
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X^Yl    ->  p{MjJd)  =  23.12  ±  (i(Ma) 

X2,  Y2    -+  p(MaUd)  =  20.70  /  fi(Ma) 

Both  the  first  and  second  principal  directions  must  be  considered  in  this  case  since 
both  represent  legitimate  solutions  to  Equation  8.5.  The  actual  vectors  from  the 
^-dimensional  subspace  that  do  provide  the  supremizing  Ud  may  be  found  by  selecting 
a  vector  Xnew  defined  as 

Xnew  =  (7  +  J8)X1  +  {a  +  jp)X2 

as  in  Equation  8.6.  Of  course  Xnew  is  a  singular  vector  of  {DMaD~l)  and  therefore 
any  scalar  multiple  of  Xnew  is  also  a  singular  vector  of  (DMaI)-1).  This  allows  the 
complex  multiplier  (7  +  j8)  of  X1  to  be  eliminated  with  no  loss  of  generality  leaving 
only  a  and  0  as  unknowns.  Vector  Xnew  and  the  corresponding  vector  Ynew  may  now 
be  written  as 

Xnew     -    X-i  +  (a  +  j/3)X2 

*  new      —  ■ 


new 


(8.7) 

and  a  2-dimensional  optimization  performed  over  a  and  /3  to  achieve  the  MPDA 
conditions  of  aligning  \Xnew\  and  \Ynew\.  Assigning  the  complex  multiplier  (a  +  j/3) 
to  X2  with  initial  conditions  of  a  =  /?  =  0  allows  this  procedure  to  be  implemented 
very  efficiently  even  when  a  /  a2.  This  assignment  serves  to  initially  weight  Xx  more 
heavily  for  cases  with  a  >  a2  while  still  providing  for  subspace  combinations  if  a  and 
a"2  coalesce. 
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Corollary  8.1:  For  systems  with  repeated  maximum  singular  value  of  multiplicity 
q  with 

inf  aiDMaD-1)  =  a2{DMaD~l)  =  ■■■=  a^DM^'1)  =  p(Ma), 
the  following  condition  must  hold: 

n(Ma)  =  p(MaUd)   |   \Xnew(at)\~\Ynew(ai)\  =  0    i  =  1,. , .  ,2(?-  1).  (8.8) 

Proof:  From  the  statement  of  the  Corollary  it  must  be  possible  to  establish  MPDA 
because  a(DMaD-x)  ^  fi(Ma)  otherwise.  This  requires  that  some  linear  combination 
of  the  first  q  columns  of  X  and  Y  respectively  must  attain  alignment.  When  alignment 
is  achieved,  the  equality  of  Equation  8.8  must  hold.  ■ 

The  resulting  g-dimensional  optimization  requires  only  vector  operations  to  align 
Xnew  and  Ynew.  Although  this  function  is  not  convex,  it  requires  only  2(q  -  1)  opti- 
mization variables  regardless  of  the  size  of  the  system  and  thus  proceeds  very  quickly. 

Example  8.2:  Using  the  system  of  Figure  and  Table  8.2,  the  solution  to  Equation 
8.8  occurs  at  a  =  -0.0009,  and  0  =  0.0921.  The  corresponding  optimal  Ud  is 

Ud  =  diag[l,  0.91  +  j0.43,  1.0  +  jO.Ol,  -.94-;0.33,  OAl+jQ.91] 

with  p(MaUd)  =  a{i)MaD-x)  =  n{Ma)  =  24.14. 

The  values  of  a  and  j3  are  not  robust  as  small  changes  in  D  cause  the  singular 
vectors  to  range  throughout  the  2-dimensional  subspace  spanned  by  Xt  and  X2.  This 
is  important  considering  that  algorithms  used  to  solve  infu  ^(DMaD'1)  depend  on  the 
analytic  derivative,  fj,  as  developed  in  Section  5.5.  Since  the  gradient  is  a  function 
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of  Xi  and  Y1,asW  approaches  a2  the  singular  vectors,  and  hence  the  gradient,  begin 
to  change  discontinuously.  This  makes  it  extremely  difficult  to  determine  whether  a 
candidate  solution  for  ft  using  only  the  gradient  at  info  a^DM^'1)  is  a  true  cusp 
or  simply  a  kiss.  Through  the  use  of  Equation  8.8  and  the  method  of  Theorem  8.1 
however,  the  system  of  Figure  8.2  is  shown  to  be  a  kiss. 

8.3     Systems  with  info -g{DMaP-1)  =  ^(DMgP-1)  ^  n(Ma) 

For  the  case  of  a  true  cusp  as  depicted  in  Figure  8.3,  it  is  not  possible  to 
rotate  \Xnew\  into  \Ynew\  because  this  condition  only  occurs  if  MPDA  can  be  attained. 
Although  no  theoretical  bound  on  the  gap  between  info  aiDMaD'1)  and  fi(Ma)  has 
been  determined,  experimental  tests  on  cusping  systems  have  revealed  disagreements 
of  more  than  14%  [43].  Determining  the  actual  value  of  ft  previously  required  solving 
for  sup^  p(MaUd)  over  n  -  1  free  variables  where  n  denotes  the  size  of  expanded 
matrix  Ma.  Since  this  computationally  intensive  optimization  is  nonconvex  the  lower 
bound  by  itself  is  not  very  useful  since  it  is  difficult  to  determine  which  of  the  local 
maxima  is  actually  fi.  Therefore,  the  upper  bound  ofaiDMal)-1)  is  generally  chosen 
as  an  estimate  for  fi(Ma)  even  though  this  leads  to  an  unnecessarily  conservative 
control  system  . 

Using  a  method  similar  to  that  of  the  previous  section,  it  may  still  be  possible 
to  find  an  Xnew  that  produces  a  Ud  (denoted  Unew)  "close"  to  the  supremizing  Ud. 
This  possibility  arises  from  the  fact  that  while  vectors  Xu...,Xq  may  individually 
change  greatly  with  small  variations  in  0MaD~l),  the  actual  subspace  spanned 
by  linear  combinations  of  X1; ,, ,  ,Xq  remains  continuous    [44,  45].    Combining  this 
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important  property  with  the  following  lemma  leads  to  an  extremely  efficient  method 
for  approximating  (i(Ma)  for  cusping  systems. 

Lemma  8.1:  For  any  A  €  CnXn  composed  of  elements  a{j,  ij  =  1, ... ,  n, 

a  lim^^A)  =  \aa\    for  fixed  i.  (8.9) 

Proof:  The  singular  values  of  A  are  the  square  roots  of  the  eigenvalues  of  AHA.  As 
an  element  along  the  diagonal  of  A  increases  or  decreases  towards  ±oo,  the  remaining 
elements  of  A  become  insignificant  in  comparison.  At  this  point  matrix  product  AH A 
may  be  approximated  by  aiiauEu  making  the  maximum  singular  value  equal  to 

Lemma  8.1  simply  shows  that  by  perturbing  a  diagonal  element  of  a  cusping 
system,  the  multiplicity  of  a  must  eventually  disappear.  In  fact,  combinations  of  the 
diagonal  elements  may  be  perturbed  providing  a  number  of  different  ways  to  break  up 
a  cusp.  This  shift  eliminates  the  cusp  and  allows  the  establishment  of  MPDA  for  the 
shifted  matrix  (denoted  Ma)  guaranteeing  that  miDa{DMaD-1)  =  /i(Ma).  MPDA 
conditions  also  provide  a  direct  method  for  computing  Ud  given  D  through  the  use  of 
Theorem  8.1.  Modifying  this  method  to  use  Xnew  and  Ynew  from  Equation  8.7  leads 
to  the  following  Proposition. 

Proposition  8.1:  Given  some  Ma  G  CnXn  for  which 

inf a(DMaD-')  =  a2{DMaD~l)  =  ■■■=  aq{DMaD~x)  ^  ^(Ma), 
the  solution  over  only  2(q  -  1)  free  variables  of 
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sup  p(MaUnew)    i  =  l,...,2(q-l)  (8.10) 

is  a  close  approximation  of  fi,(Ma). 

Discussion:  The  continuity  of  the  spanning  subspace  Xt  •  •  •  Xq  implies  that  small 
perturbations  in  the  matrix  elements  of  Ma  will  result  in  correspondingly  small  pertur- 
bations in  the  subspace.  As  a  scalar  indication  of  the  distance  between  two  matrices 
A, B  €  Cnxn,  define  an  error  function  E  =  \\A  -  B\\2  [42].  Shifting  the  diagonal 
elements  of  the  cusping  system  in  Table  8.3  until  the  cusp  disappears  produces  a  new 
system  with  desirable  MPDA  properties.  If  this  new  system  is  sufficiently  close  to 
the  original  cusping  system  then  Xnew  of  the  two  systems  should  be  similar. 

Examination  of  Tables  8.2  and  8.3  reveals  that  the  two  systems  are  identical 
except  for  an  additive  shift  of  +14  to  all  the  diagonal  elements  of  Table  8.2  so  that 
Ma  =  Ma  +  14  ■  h.  Denote  Ma  and  the  infimizing  D  of  Table  8.2  as  Ma2  and  D2 
respectively.  Similarly  denote  the  corresponding  values  of  Table  8.3  as  Ma3  and  D3. 
Then 


E  =  \\D3Ma3D^  -  D2Ma2D21\\  =  16.44. 

Since  a(D3Ma3D^)  =  13.11  from  Figure  8.3,  E  indicates  that  the  two  systems  are 
not  particularly  close. 

However,  an  alternative  shift  of  only  element  M„(l,  1)  by  +4  also  breaks  up  the 
cusp.  Figure  8.4  and  Table  8.4  show  this  new  system  with 

H{Ma)  =  a{bAMaiD4l)  =  15.182. 
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13.5 


1.5       1.6       1.7       1.8       1.9         2         2.1       2.2       2.3       2.4       2.5 

D(l,l) 


Figure  8.4:       System  from  Figure  8.3  with  Shifted  Af„(l,  1)  Element. 


9.18  +  0.37; 
-0.20  -  3.07; 

6.42  +  1.85; 

0.06  +  0.64; 
-2.39  -  5.34; 


Table  8.4:  Ma  Matrix  for  Figure  8.4. 


6.82  -  1.75; 

4.56  +  1.29; 
-0.70  +  1.03; 
-0.53  -  2.47; 

3.21  -  0.78; 


3.13  -  0.95; 

-1.44  +  0.35; 

1.34-2.01; 

3.53  -  0.97; 

3.74  +  1.38; 


-4.92  +  1.11; 

3.22  +  2.37; 
-0.77  -  0.82; 
-3.03  -  3.11; 

3.24  -  0.03; 


3.34  -  4.59; 
-1.32  +  3.15; 
-0.13  +  1.36; 

3.93  -  0.96; 

1.41-0.33; 
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Error  function  E  then  becomes 


E  =  \\D3Ma3D^  -  btM^blH  =  4.77 


indicating  that  System  4  is  much  closer  to  System  3  than  System  2  was.  Since 
the  spanning  subspace  is  continuous,  it  seems  reasonable  for  the  optimum  Xmw  = 
Xx  +  («i  +  ja2)X2  of  the  cusping  system  to  be  similar  to  that  of  Table  8.4  for  which 
MPDA  is  established.  Table  8.5  shows  the  variation  of  Xnew  as  the  cusping  system 
of  Table  8.3  (Shift=0)  is  shifted  to  the  noncusping  system  of  Table  8.4  (Shift=4). 

Table  8.5:  Variation  of  Xnew  due  to  Shift  of  Table  8.3  System. 

Shift  of  Ma(l,l) 


+0 


-0.053-  j  0.237 
0.425+  j  0.210 
0.300+  j  0.240 
-0.027-  j  0.345 


+1 


-0.070  -  j  0.248 
0.424  +  j  0.166 
0.309  +  j  0.246 
-0.058  -  j  0.352 


+2 


-0.080  -  j  0.254 
0.409  +  j  0.135 
0.314  +  j  0.244 
-0.075  -  j  0.355 


+3 


-0.089  -  j  0.258 
0.388  +  j  0.112 
0.315  +  j  0.241 
-0.087  -  j  0.354 


+4 


1 


-0.10  -  j  0.26 
0.37  +  j  0.09 
0.31  +  j  0.24 
-0.09  -  j  0.35 


Since  shifting  other  combinations  of  matrix  elements  will  also  break  up  the  cusp, 
it  may  certainly  be  possible  to  find  noncusping  systems  for  which  E  is  less  than  the 
4.77  achieved  by  System  4.  Regardless  of  the  shifting  combination  used  to  break  up 
a  cusp,  the  2(q  -  l)-dimensional  subspace  spanned  by  Xx,...  ,Xq  must  vary  contin- 
uously from  the  cusping  condition  to  any  of  the  noncusping  conditions.  This  implies 
that  the  spanning  subspace  is  sufficiently  well-behaved  to  provide  an  Xnew  from  Equa- 
tion 8.7  and  corresponding  Unew  that  closely  approximates  Ud  of  p(MaUd)  =  /i(Ma). 
Although  no  bounds  have  been  placed  on  the  difference  between  supa[  p(MaUnew)  and 
sup^  p(MaUd),  numerical  experience  suggests  that  this  difference  is  quite  small.   In 
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fact,  the  largest  deviation  between  this  estimate  and  fi  found  to  date  is  less  than 
0.34%.  ■ 

Applying  Equation  8.10  to  the  cusping  system  of  Table  8.3  gives 

Unew  =  diag[l,  0.75  -  j0.67,  -.48  +  j0.88,  -1.00-j0.02,  0.78-j0.63] 

with  «!  =  -.23,  a2  =  -0.11  and  p{MaUnew)  =  12.768. 

To  determine  whether  Unew  is  actually  a  candidate  for  //,  the  derivative  dp/d&i  for 
l  =  1, . . , ,  n  at  Unew  may  be  examined  for  its  proximity  to  the  0  vector.  In  this  case 
it  is  found  to  be 

0.24 


dp_ 


0.10 

-0.19 

0.31 

0.17 


indicating  that  Unew  is  indeed  "almost"  a  stationary  point.  Using  Unew  as  an  initial 
condition,  a  direct  supremization  over  all  Ud  returns  a  value  of 

Ud  =  diag[l,  0.59-j0.80,  -.53  +  ;0.85,  -.99-jO.ll,  0.66-j0.75]  (8.11) 


with  p(MaUd)  =  n{Ma)  =  12.81.  The  difference  between  p(MaUnew)  and  p(MaUd)  in 
this  case  is  0.33%. 

This  particular  example  was  chosen  because  it  produced  the  largest  percent  dif- 
ference between  p{MaUnew)  and  p,(Ma)  out  of  over  100  cusping  systems  tested.  The 
next  largest  percent  difference  was  10  times  smaller  than  this  case:    about  0.03%. 
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While  larger  differences  are  certainly  possible,  the  fact  that  in  the  majority  of  ex- 
amples investigated,  p{MaUnew)  and  p(MaUd)  agreed  within  four  significant  figures 
is  most  encouraging.  Table  8.6  on  page  114  contains  a  comparison  between  fi(Ma) 
and  p(MaUnew)  for  several  cusping  systems.  (This  comparison  will  be  discussed  more 
completely  in  Section  8.8.) 

8.4     Effects  of  Nonconvexity 

It  must  be  emphasized  that  both  the  2(q  -  1)- dimensional  search  over  the  a,- 
and  the  n  -  1  dimensional  search  over  Ud  are  nonconvex.  Algorithms  designed  to 
find  the  global  maximum  must  employ  some  procedure  to  search  the  domains  of  Ud 
or  Unew  respectively.  Therefore,  the  advantages  of  the  2 (q  -  1)- dimensional  search 
become  more  profound  because  the  domain  of  Unew  does  not  increase  with  system 
size. 

Although  repeated  singular  values  with  multiplicity  q  =  n  are  theoretically  possi- 
ble, by  far  the  more  common  case  is  for  q  =  2.  For  such  systems,  an  additional  bonus 
of  a  2(q  -  1)  =  2-dimensional  search  is  that  3-dimensional  mesh  plots  of  p(MaUd) 
vs  qj  and  a2  may  be  generated  allowing  all  maxima  to  be  examined.  This  effectively 
eliminates  the  difficulties  of  multiple  maxima  because  all  candidate  points  for  p  can 
easily  be  examined. 

Such  a  mesh  plot  for  the  system  of  Figure  8.3  is  presented  in  Figure  8.5.  In 
this  plot,  at  and  a2  range  from  -2  to  2.  Two  pronounced  local  maxima  are  readily 
apparent:  one  corresponding  to  the  value  of  p(MaUnew)  =  12.768  and  the  other  at 
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<l.^  <l^ 


Figure  8.5:       Mesh  Plot  of  System  from  Table  8.3. 
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p{MaUnew)  =  12.741  where 

Unew  =  diag[l,  -.13  +  j0.99,  1.00  +  j0.084,  -.91+j0.42,  -0.81+^0.58] 

with  ax  =  0.25  and  a2  =  0.54.   Using  Unew  as  a  starting  point,  a  search  over  all  Ud 
results  in  a  local  stationary  point  at 

Ud  =  diag[l,  -.17  +  ^0.98,  0.95+j0.32,  -.94  +  ^0.33,  -.92+^0.40] 
with  corresponding  p(MaUd)  =  12.76  ^  ft(M9)  =  12.81. 

8.5     Range  of  o?t- 

The  ability  to  easily  select  local  maxima  of  p(MaUnew)  from  a  mesh  plot  as 
candidates  for  fi  is  a  tremendous  improvement  over  previous  methods  which  require 
optimizing  p(MaUd)  over  as  many  as  n  -  1  variables,  especially  as  the  system  size 
increases.  The  range  of  <*i  and  a2  used  in  Figure  8.4  appears  to  adequately  cover 
the  subspace  spanned  by  Xx  and  X2  allowing  Unew  to  be  found.  As  ax  and  a2 
move  away  from  the  local  maxima,  the  magnitude  of  p  decreases  suggesting,  (but  not 
guaranteeing),  that  no  additional  points  of  interest  exist  outside  of  the  plotted  region. 
It  is  important  to  determine  the  maximum  range  of  the  a,-  because  if  the  bounds  can 
only  be  placed  at  ±oo,  then  a  mesh  plot  or  any  other  search  method  would  clearly 
be  impractical. 

For  the  case  of  multiplicity  q  =  2,  determining  limits  on  ax  and  a2  requires  ex- 
amination of  Equation  8.7.  Placing  the  complex  multiplier  (ctj  +  ja2)  on  X2  makes 
the  worst  case  condition  when  the  optimum  Xnew  =  X2  so  that  Xx  makes  no  con- 
tribution.  Such  a  condition  implies  that  ai  =  oo  is  required  to  produce  the  desired 
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Xnew  —  X2.  However,  since  X\  and  X2  are  singular  vectors  they  may  be  normalized 
such  that  \\Xi\\  =  \\X2\\  =  1.  Since  the  2— norm  of  a  vector  is  greater  than  or  equal  to 
any  element  of  that  vector,  no  element  of  X\  or  X2  can  be  greater  than  1.  Therefore, 
in  forming  Xnew  =  X\  +  {ot\  +  ja2)X2,  a.\  and/or  a2  values  much  greater  than  1  will 
essentially  eliminate  any  contribution  of  X\  so  Xnew  quickly  reaches  its  asymptotic 
value  of  X2.  In  the  rare  case  that  such  a  worst  case  condition  occurs,  values  of  a^  or 
a2  of  20  or  less  should  be  more  than  adequate  to  produce  the  optimum  Unew.  In  prac- 
tice, large  values  of  cki  and/or  a2  can  be  eliminated  by  merely  shifting  the  complex 
multiplier  from  X2  to  X\  so  that 

Xnew  =  («i  +  ja2)X1  +  X2 

For  Xnew  ~  X2,  the  condition  of  c*i  =  a2  =  0  then  prevents  any  contribution  from 
Xi. 

8.6     Principal  Direction  Alignment 

The  ability  to  compute  candidates  for  fi  using  only  2(q  —  1)  variables  provides  an 
extremely  practical  means  for  dealing  with  cusping  singular  values.  After  determining 
supa.  p(MaUnew)  the  result  may  compared  to  the  singular  value  upper  bound  and  if 
the  conservatism  gap  does  not  exceed  desired  tolerances  no  further  optimization  is 
required.  The  conservative  upper  bound  is  simply  chosen  as  an  estimate  for  fi  and 
the  frequency  sweep  over  the  desired  range  continues.  In  the  event  that  gap  between 
the  upper  and  lower  bounds  exceeds  the  termination  requirements,  the  supremizing 
matrix  Unew  offers  an  improved  initial  condition  for  a  search  over  the  full  set  (n  —  1)  of 
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optimization  variables.    Unfortunately,  the  computational  aspects  of  this  additional 

search  are  potentially  overwhelming:  especially  as  the  system  size  increases.  Rather 
than  continuing  the  search  for  fi  using  eigenvalue  calculations,  it  will  be  shown  that  a 
simple  transformation  of  Unew  allows  the  optimization  to  proceed  using  singular  value 
decompositions.  Because  the  singular  values  are  the  square  roots  of  the  eigenvalues 
of  positive  definite,  symmetric  matrices,  (AHA  or  AAH),  the  singular  value  decompo- 
sition offers  significant  improvements  in  computing  efficiency  over  general  eigenvalue 
decompositions.  For  example,  using  the  MATLAB™  "eig"  and  "svd"  functions  on 
10  random,  complex  9x9  matrices,  the  average  flop  counts  are  as  follows: 

eig  — >  53,606  flops 

svd^  38,448  flops. 

As  the  system  size  increases,  this  difference  becomes  even  more  significant  indicat- 
ing the  obvious  computational  advantage  of  the  singular  value  decomposition  versus 
the  eigenvalue  decomposition.  It  must  also  be  noted  that  the  flop  counts  above  in- 
clude calculation  of  the  singular  vectors  while  the  eigenvectors  were  not  similarly 
determined. 

The  MPDA  property  of  Chapter  4  shows  that  the  unique  stationary  point  as- 
sociated with  info  a{DMaD~l)  corresponds  to  a  stationary  point  of  supjy  p(MaUd)- 
When  a  repeats  there  may  be  no  stationary  point  of  ~&  so  the  direct  relationship  be- 
tween a  and  [i  is  in  general  lost  (the  "kiss"  condition  is  an  exception).  However,  using 
the  principal  direction  alignment  (PDA)  property  developed  by  Daniel,  Kouvaritakis 
and  Latchman  ,  a  similar  relationship  is  shown  to  exist  between  Gi{DMaD~x)  and 
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p(MaUd)  where  the  <j;  represents  some  singular  value  other  than  the  maximum  [31]. 
Development  of  the  PDA  property  starts  with  the  following  lemma  describing  the 
well  known  eigenvalue  shift  property. 

Lemma  8.2:  For  any  matrix  A  €  CnXn  and  scalar  r 

Xi(A  +  rI)  =  Xi(A)  +  r 

Proof:  Denote  the  eigenvalues  of  matrix  A  as  A  and  of  matrix  (A  +  rl)  as  A.  The 
respective  eigenvalues  are  simply  the  solutions  to 

det{XI-A}  =  0 

and 

det{XI  -  A  -  rl}  =  0. 
Now  define  (A  —  r)  =  A  and  rewrite  the  second  equation  as 

det{XI~  A]  =  0. 

Since  the  eigenvalues  of  a  matrix  are  the  roots  of  its  characteristic  equation  they 
must  be  unique;  therefore  A  =  A.  Substituting  back  for  A  gives  A  =  A  +  r  and  thus 
completes  the  proof.  ■ 

While  the  eigenvalues  of  a  matrix  obey  the  shift  property,  the  singular  values  do 
not.  As  in  Section  8.3,  this  provides  a  means  of  breaking  up  the  repeated  singular  val- 
ues by  introducing  a  sufficiently  large  r.  The  spectral  radius /singular  value  inequality 
then  becomes 

p(MaUd)  +  r  =  p(MaUd  +  rl)  <  mfaiDMMdD-1  +  rl).  (8.12) 
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Equation  8.12  requires  that  Xz{MaUd)  corresponding  to  p(MaUd)  be  real.    Such  a 

requirement  is  not  a  limitation,  however,  since  Ud  may  be  rotated  by  a  scalar  multiplier 

e*  without  affecting  the  magnitude  of  p{MaUd).    Rewriting  Equation  8.12  into  a 

structured  singular  value  relationship  results  in 

fi(Ma)  =  sup  p(MaUd)  =  sup  p(MaUd  +  r/)  -  r  =  sup  mi  afDMa^D'1  +  rl)~r 

ud  Ud  Ud     D 

(8.13) 
Unfortunately,  the  convexity  properties  that  made  the  optimization  infD  ^{DMaD'1) 
so  attractive  originally  are  now  lost  because  the  right-hand  side  of  Equation  8.13  is 
no  longer  invariant  to  Ud.  Still,  the  addition  of  scalar  r  offers  an  alternative  means 
of  establishing  MPDA  for  cusping  systems.  The  following  theorem  is  originally  due 
to  Latchman    [27]. 

Theorem  8.2:  Let  a(DMaUdD-1  +  rl)  have  stationary  points  at  Ud  and  D  respec- 
tively. If  r  is  large  enough  to  ensure  that  cr^M^I)-1  +  rl)  is  simple,  then  MPDA 
is  established  for  the  matrix  (DMaUdD-1  +  rl)j  where  J  is  a  matrix  containing  only 
±1  along  the  diagonal  and  zeros  elsewhere. 

Proof:  The  stationary  points  of  sup^inf^  a{DMa\JdD~x  +  rl)  occur  when 

^[5(Map-'  +  r/)]  =  0  (8.14) 


anc 


^\a{DMaUdD~l  +  rl)}  =  0.  (8.15) 

The  stationary  points  of  Equation  8.14  have  already  been  shown  to  occur  upon  es- 
tablishment of  MPDA  for  the  matrix  {DMaUdD~l  +  rl).  This  condition  applies  only 
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when 

MQi|  =  lKi|  (8.16) 

for  all  i.  Differentiation  of  Equation  8.15  follows  in  a  manner  similar  to  that  of  Section 
5.5  so  that  the  stationary  points  are  characterized  by 

Jj{Y1HEiiXi-X*IEiiY1]  =  0 

or 

Im{X?EuYi)  =  0  (8.17) 

for  all  i.  The  requirement  of  Equation  8.17  may  be  stated  in  terms  of  the  arguments 
of  the  individual  vector  elements  such  that 

axg(X.-i)  =  arg(Ki)  +  mi7r  (8.18) 

for  all  i  with  integer  m,-.  Applying  Equations  8.16  and  8.18  requires  that 

X,  =  JYX 

and  completes  the  proof.  ■ 

Progressing  from  this  point  to  MPDA  therefore  requires  that  J  =  /  for  all  cases  in 
which  both  Equations  8.14  and  8.15  are  satisfied.  The  following  theorem,  originally 
due  to  Latchman  with  simplifications  by  Young    [27,  46]  shows  that  J  must  equal  /. 

Theorem  8.3:  Let  ^{DMJJ&D^1  +  rl)  have  stationary  points  at  Uj.  and  D  respec- 
tively. If  r  is  large  enough  to  ensure  that  a(DMaUdb~1  +  rl)  is  simple,  then  MPDA 
is  established  for  the  matrix  (DMaUjD"1  +  rl). 
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Proof:  As  mentioned  above,  this  proof  reduces  to  showing  that  J  =   I  at  all 
stationary  points  of  Equation  8.13.   Matrix  {DMaUdD'1  +  rl)j  has  been  shown  to 
satisfy  the  MPDA  requirement.  Therefore 

p{DMaUdI)-1  +  rl)  <  a{\DMaUdb-1  +  rl]j). 
Multiplying  by  Xf  and  X%  gives 

X^DMaUdD-1  +  rI)X1  <  X?{DMaUdb-'  +  rI)JXt 
which  may  be  simplified  as 

rX? (I  -  J)X1  <  X?DMab-lUd{J  -  I)X1.  (8.19) 

Clearly  as  r  increases,  the  left-hand  side  of  Equation  8.19  increases  while  the  right- 
hand  side  remains  unchanged.  The  inequality  would  therefore  be  contradicted  unless 
J  =  I  which  completes  the  proof.  ■ 

Theorem  8.2  reveals  that  for  r  large  enough  to  separate  a  and  a2 

svLj>p(MaUd)  =  supmiaiDMa^D-1  +  rl)  -  r. 
ud  ud    d 

More  importantly,  the  following  lemma  serves  to  relate  stationary  points  of 
suP[/dinfD  v{DMaUdD~l  +  rl)  to  those  of  ^MaUdD-1)  where  the  trv,  represent 
singular  values  with  magnitudes  less  than  a.  This  new  formulation  is  invariant  to  Ud 
and  no  longer  requires  the  shift  by  r. 

Lemma  8.3:  For  D  and  Ud  defined  as  optimizing  solutions  to  Equation  8.13,  denote 
Xx  as  the  left  major  singular  vector  of  {bM$J)-*+rT).  Then  Xx  and  its  correspond- 
ing conjugate  transpose,  Xf ,  are  also  right  and  left  eigenvectors  of  (DMaUdb-1).  In 


-1 

r          -, 

■ 

Ai 

0 

x? 

Xi 

w 

0 

A 

yT 
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addition,  matrix  {bM^db'1)  has  some  right  singular  vector  which  is  aligned  with 
its  corresponding  left  singular  vector.  This  alignment  between  nonmajor  singular 
vectors  will  be  referred  to  as  Principal  Direction  Alignment    [31]. 

Proof:  From  Theorem  8.2  it  is  known  that  (DM^D-1  +  rl)  has  MPDA  and 
therefore  Xr  and  X?  form  a  left  and  right  eigenvector  pair  for  (DMaUdb'1  +  rl), 
From  the  eigenvalue  shift  property  of  Lemma  8.2,  Xx  and  Xf*  must  form  a  left  and 
right  eigenvector  pair  for  (DMaUdf)-1)  as  well.  This  allows  the  eigenvalue  decompo- 
sition of  (DMaUji)-1)  to  be  written  as 


DMMdD-1  = 


where  A,  W  and  VT  include  all  the  eigenvalues,  and  right  and  left  eigenvectors  except 
Ai ,  X\  and  Xf  respectively.  Inspection  of  the  hermetian  products 

{DMaUib-l)H{DMaUdb-1) 

and 

{bMaudb-l){bMaudb-l)H 

as  in  Equation  8.4,  reveals  that  Xi  is  both  a  right  and  corresponding  left  singu- 
lar vector  associated  with  some  a^DMaUdb-1).  Therefore,  (DMaUdb-1)  has  PDA 
concluding  the  proof.  ■ 

Combining  the  results  from  Theorem  8.2  and  Lemma  8.3  leads  to  the  following 
theorem  concerning  the  direct  relationship  between  sup^  p{MaUd)  and 
^(DMaUdb-1). 
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Theorem  8.4:  For  any  matrix  Ma,  there  always  exists  an  optimal  U4  and  D  such 

that  (DMaUdD'1)  achieves  PDA  and  the  corresponding  singular  value  is  a  stationary 

point  with  respect  to  D.    Also,  under  these  circumstances  the  following  condition 


arises: 


sup  p{MaUd)  =  maxj^pMaZr1)  |  ~[ai(DMaUdD-1)]  =  ol  (8.20) 

ud  *      {  odj  J 

for  i,j  =  l,...,n. 

Proof:  Theorem  8.2  guarantees  the  existence  of  an  optimizing  Ud  and  D  such  that 
(DMaUdD'1)  satisfies  MPDA.  The  shift  property  may  then  be  employed  giving 

p{MaUd)  =  a{DMaUdb-1  +rl)-r  =  ^{DM^-1). 

Theorem  8.3  shows  that  for  optimal  D  and  Ud  matrix  (DMaUdb-1)  will  have  PDA. 
For  the  case  of  i  =  1,  a  does  not  cusp  and  the  MPDA  property  is  achieved.  When 
i  >  1,  achieving  PDA  for  (DMaUj,!)-1)  requires  one  of  the  singular  vector  pairs  to  be 
aligned  indicating  that  ai(DMaD~l)  is  a  stationary  point  with  respect  to  D.  Selecting 
the  largest  such  stationary  point  completes  the  proof.  ■ 

Example  8.3:  This  example  illustrates  these  results  using  the  system  of  Table  8.3. 
The  supremizing  Ud  is  the  same  as  that  presented  in  Equation  8.11.  Infimizing  the 
right-hand  side  of  Equation  8.20  over  D  with  r  >  1.87  gives 

D  =  diag[l,  0.749,  1.106,  1.225,  1.218] 


and 


fi{Ma)  =  p{MaUd)  =  a{DMaUdb-1  +rl)-r  =  a2{DMaD-1)  =  12.81. 
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The  relationship  of  Equation  8.20  allows  the  singular  value  decomposition  to  be 

retained  in  the  computation  of  ft  even  for  repeated  a.  Unfortunately  the  stationarity 
condition  on  a{  no  longer  restricts  the  search  for  p  to  minima  or  maxima  but  to  inflec- 
tion points  as  well.  This  introduces  additional  complexity  to  an  already  nonconvex 
function  as  the  search  criteria  now  includes  sorting  through  all  stationary  points.  If, 
however,  the  search  over  the  a;  produces  a  value  of  p(MaUnew)  close  to  p(Ma),  then 
an  optimization  based  on  the  right-hand  side  of  Equation  8.20  may  offer  tremendous 
computational  advantages  over  the  left-hand  side  of  the  equation. 

8.7     Direct  Calculation  of  D  from  Ud 

Section  8.6  introduced  a  method  for  directly  calculating  Ud  given  D  such  that 
aiDMal)-1)  =  p(MaUd)  =  p(Ma).  Using  the  2(q  -  l)-dimensional  search  over  az, 
a  candidate  for  Ud  denoted  Unew  may  be  calculated.  In  order  for  this  estimate  to 
serve  as  a  useful  initial  guess  for  the  right-hand  side  of  Equation  8.20,  some  means  of 
transforming  Unew  into  a  corresponding  Dnew  must  be  found.  Such  a  transformation 
was  originally  presented  by  Fan  and  Tits  as  a  means  of  verifying  that  the  solution 
to  their  nonconvex  vector  optimization  method  equals  p  [38].  By  modifying  their 
original  theorem  and  proof  pertaining  to  this  transformation  it  is  possible  to  directly 
determine  D  corresponding  to  Ud  such  that  a^DM^-^)  =  p{MaUd)  even  when  a 
repeats.  Equation  7.2  is  repeated  here  for  convenience. 

tx{Ma)  =  max{||Afaar||a  |  ||P,-a:||2||M0x||2  =  \\VtMax\\2,  i  =  l,...,m}  (8.21) 


Ill 

The  following  proposition  and  subsequent  theorem  are  originally  due  to  Fan  and 
Tits  with  modifications  to  account  for  systems  with  repeated  singular  values    [38]. 

Proposition  8.2:  Suppose  x  6  Cn  is  a  solution  to  Equation  8.21.  Then  there  exists 
a  vector  A;,    i  =  1, . . . ,  m  such  that  the  following  equation  holds: 

m 

M?Max  +  £  A,-  [M»ViMax  -  \\V%x\\2M?  Max  -  ||Max||22^x]  =  0  (8.22) 

8  =  1 

where  the  Vz  are  projection  matrices  as  in  Chapter  7  and  the  A\-,s  are  multipliers  with 
A;  >  0,  i  =  1,. . .  ,  n. 

Theorem  8.5:  Let  A  =    Ai , . . . ,  A„j  be  the  unique  multiplier  vector  corresponding 


to  the  optimizing  x.  Then  D  —  diag 


i 


\  2  \   ■ 

A\     , .  . .  ,  An 


is  the  solution  to  Equation  8.20. 


Proof:  From  Chapter  7,  some  Y{  must  exist  such  that 

b-1M^D2MaD-lYi  =  aiiDMaD-^Yi  =  »2(Ma)Yi  =  \\Max\\22Yt. 
Multiplying  on  the  left  by  D  gives  the  expression 

M»D2Max  =  \\Max\\22D2x.  (8.23) 

Now  define  some  f3  such  that 

* 2 
(3i  =  adi     i  =  l,...,m  (8.24) 

with  a  selected  such  that  ££,  Pi\\Vix\$  =  1.   Then  at)2  =  ££,  fcVi  and  Equation 
8.23  becomes 

E  fii  \M?ViMax  -  \\Max\\lVix\  =  0 


giving 


n 

M?Max  +  J2&  [M?ViMax  -  \\ViX\\lM?Max  -  \\Max\\22Vtx}  =  0. 


i=i 
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Equating  $  =  A;,  this  equation  becomes  identical  to  Equation  8.22  and  through  the 
relationship  of  Equation  8.24  the  proof  is  complete.  ■ 

It  is  now  possible  to  directly  determine  a  D  given  Ud  such  that 
p(MaUd)  =  ai(DMab-1)  =  fi{Ma).   Again  using  Ud  from  Equation  8.11,  D  is  found 
to  be 

D  =  diag[l,  0.749,  1.106,  1.225,  1.218] 

with 

fi(Ma)  =  p(MaUd)  =  <j2(DMab-1)  =  12.81 

providing  the  same  results  as  before.  However,  no  additional  optimization  is  required 
for  this  method  as  opposed  to  Example  8.3  which  requires  an  infimization  over  D  of 
^{DMaUciD'1  +  rl).  Provided  that  the  2(q  -  l)-dimensional  optimization  generates 
a  Unew  close  to  Ud,  the  transformation  from  Unew  to  Dnew  may  be  invoked  to  solve 
for  p  using  the  method  of  Equation  8.20. 

8.8     Comparison  of  Optimization  Methods 

The  results  of  this  chapter  are  summarized  in  Tables  8.6  and  8.7.  Table  8.6 
lists  system  size,  the  results  of  ixdDu{DMaD-1),  the  value  fl(Ma),  the  results  of 
SUP«, ;p(MaUnew),  and  the  per  cent  difference  between  fi(Mu)  and  supa.  p{MaUnew). 
System  6  corresponds  to  that  of  Figure  8.3  and,  as  mentioned  earlier,  it  shows  the 
largest  deviation  between  (i(Ma)  and  supat  p(MaUnew)  of  all  cusping  systems  studied 
to  date.  System  5  shows  the  largest  deviation  between  miDa(DMaD~1)  and  p,{Ma) 
while  the  corresponding  value  of  supQ.  p{MaUnew)  equals  p,{Ma)  within  four  significant 
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figures.  Clearly  the  result  of  the  2(q  -  1)- dimensional  optimization  supa.  p{MaUnew) 
provides  an  excellent  estimate  of  \i. 

Depending  on  the  particular  control  system  under  development,  it  is  likely  that 
these  estimates  of  (j,  will  be  satisfactory  without  any  additional  optimization.  In  the 
event  that  ft  must  be  determined  exactly,  the  optimum  Unew  from  sup„.  p(MaUnew) 
may  be  used  as  an  starting  point  for  further  searches.  Using  the  same  systems  as 
Table  8.6,  Table  8.7  summarizes  the  flop  count  requirements  of  the  various  methods. 
Each  of  the  optimization  methods  call  a  series  of  MATLAB™  m-files  developed 
by  Grace  [33].  While  m-files  generally  execute  more  slowly  than  compiled  code,  the 
ability  to  track  flop  counts  provides  a  means  of  comparing  the  efficiency  of  algorithms 
regardless  of  the  host  machine. 

In  order  to  provide  as  equitable  a  comparison  between  methods  as  possible  each 
cusping  system  was  required  to  satisfy  three  criteria.  Of  all  cusping  systems  tested, 
only  those  meeting  these  criteria  are  included  in  Tables  8.6  and  8.7.  (Note:  It  should 
be  mentioned  that  failure  to  satisfy  these  criteria  does  not  mean  that  the  optimization 
methods  could  not  find  (i  for  a  particular  system,  it  simply  means  that  any  flop  count 
comparison  would  be  erroneous.) 

The  first  selection  criterion  is  that  optimization  methods  A  and  B  (sup^  p(MaUd) 
and  sup„.  p(MaUnew)  respectively),  must  begin  at  the  same  starting  point.  This  point 
is  a  value  of  Ud  generated  by  the  method  of  Theorem  8.1.  For  q  =  2,  two  values  of 
Ud  are  generated  with  Ux  corresponding  to  Xx  and  Yi  and  U2  corresponding  to  X2 
and  Y2.   Since  these  singular  vector  pairs  are  effectively  interchangeable  at  the  cusp 
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Table  8.6:  Cusping  Systems 


%  Difference 

System 

Size 

4 

infv  a(DMaD-1) 
6.8125 

6.7803 

supa,  p{MaUnew) 

(i  vs  p(MaUnew) 

1 

6.7803 

-0.000 

2 

4 

5.9398 

5.9015 

5.9015 

-0.000 

3 

4 

7.1238 

6.9697 

6.9488 

-0.013 

4 

4 

6.5750 

6.3673 

6.3671 

-0.003 

5 

4 

1.0322 

0.9470 

0.9470 

-0.000 

6 

5 

13.114 

12.810 

12.761 

-0.326 

7 

5 

31.710 

31.443 

31.440 

-0.009 

8 

5 

35.838 

35.615 

35.604 

-0.031 

9 

5 

43.913 

43.809 

43.807 

-0.007 

10 

5 

35.126 

34.930 

34.922 

-0.022 

11 

9 

7.7400 

7.7394 

7.7394 

-0.000 

12 

9 

7.8366 

7.8346 

7.8346 

-0.000 

13 

9 

9.7072 

9.6193 

9.6185 

-0.008 

14 

9 

29.473 

29.254 

29.253 

-0.003 

15 

9 

128.02 

127.07 

127.07 

-0.002 

16 

9 

74.749 

74.502 

74.501 

-0.001 

17 

9 

151.95 

150.29 

150.28 

-0.002 

18 

9 

209.64 

209.21 

209.20 

-0.001 

19 

9 

35.238 

34.610 

34.610 

-0.001 

20 

9 

35.918 

35.339 

35.336 

-0.008 
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Table  8.7:  Floating  Point  Operations  Comparison 


Optimization  Method 

A 

B 

C 

D 

B  [ 

B  i 

System 

supu,  p(MaUd) 

supa.  p(MaUnew) 

suPt/w  p(MaUd) 

m&XiCr^DMaD-1) 

1 

234,427 

128,578 

149,067 

149,899 

2 

138,177 

113,174 

144,883 

150,354 

3 

128,398 

136,015 

191,239 

157,372 

4 

188,109 

193,382 

212,596 

214,373 

5 

171,721 

63,162 

74,317 

79,787 

Average: 

172,166 

126,862 

154,420 

150,357 

6 

579,929 

301,294 

395,592 

343,173 

7 

382,910 

299,005 

550,583 

343,313 

8 

709,431 

285,912 

452,255 

329,634 

9 

839,772 

553,011 

722,400 

595,733 

10 

415,088 

256,610 

582,077 

300,868 

Average: 

585,426 

339,166 

540,581 

382,544 

11 

677,513 

226,805 

408,046 

452,390 

12 

1,037,283 

384,400 

563,199 

611,346 

13 

1,992,286 

816,698 

1,221,836 

1,039,923 

14 

2,909,386 

1,172,394 

1,850,040 

1,398,542 

15 

2,850,485 

1,078,454 

2,345,768 

1,307,238 

16 

3,203,488 

984,491 

1,952,090 

1,211,701 

17 

3,327,218 

1,240,607 

2,200,240 

1,465,977 

18 

2,786,795 

1,330,380 

2,342,328 

1,556,346 

19 

3,349,508 

946,319 

1,687,675 

1,113,345 

20 

4,492,976 

1,684,968 

3,521,077 

1,937,664 

Average: 

2,662,694 

986,551 

1,809,229 

1,109,472 
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(see  Equations  8.4  and  8.5),  both  Ui  and  U2  provide  candidate  starting  points  for 
the  subsequent  optimizations.  The  choice  between  the  two  is  made  by  selecting  the 
maximizer  of  p(MaUi)  for  %  =  1, . . , ,  q,  The  extra  flops  required  to  select  the  desired 
Ui  are  generally  more  than  offset  by  increased  convergence  rates  resulting  from  the 
superior  starting  point. 

Choosing  U  or  U2  as  a  starting  point  for  sup^  p(MaUd)  corresponds  to  choosing 
whether  the  complex  multiplier  (aj  +  ja2)  acts  on  Xx  or  X2  in  Equation  8.7.  Using 
<*i  =  ol2  =  0  as  an  initial  condition,  matrix  U\  results  from 

Xnew  —  X\  +  (ax  +  ja2)X2 
while  U2  results  from 

Xnew  =  («1  +  JOi2)Xi  +  X2. 

This  procedure  allows  sup^  p{MaUd)  and  sup^  P{MaUnew)  to  start  at  the  same  point 
satisfying  the  first  selection  criterion. 

The  second  criterion  requires  the  optimization  to  converge  to  fi  from  the  common 
starting  point.  This  greatly  reduces  the  number  of  systems  available  for  flop  count 
comparisons  because  the  methods  are  all  known  to  be  nonconvex.  While  changing 
the  initial  starting  points  would  eventually  lead  to  //,  the  possibility  of  an  equitable 
flop  count  comparison  is  eliminated.  This  local  convexity  criterion  extends  to  the 
methods  of  columns  C  and  D  of  Table  8.7  as  well.  For  these  methods,  supQt  p(MaUnew) 
provides  an  improved  starting  point  for  subsequent  optimizations  by  sup^  p(MaUd) 
and  m&Xiai{DMaD-1)  respectively.  As  discussed  in  Section  8.6,  the  optimization 
based  on  ^(DMaD'1)  searches  for  all  stationary  points  including  maxima,  minima, 


1  <  0.02.    Corresponding  parameters  for  max;  <?i(DMaD~x)  were  then  set  such 
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and  inflection  points.  While  the  computational  advantages  of  the  singular  value 
decomposition  makes  it  quite  fast,  its  non convexity  generally  prevents  its  use  without 
the  improved  starting  point  provided  by  supa.  p(MaUneU!). 

The  final  selection  criterion  involves  termination  conditions  for  the  separate  meth- 
ods. While  the  end  results  of  methods  A,  C,  and  D  should  be  identical,  differences  in 
problem  formulation  prevent  the  use  of  a  common  termination  condition.  If  one  algo- 
rithm employs  a  more  strict  termination  condition,  it  may  unnecessarily  require  more 
flops  than  a  second,  less  strict  method  to  achieve  the  same  results.  To  resolve  this 
dilemma,  parameters  were  adjusted  to  halt  the  optimization  of  sup^  p(MaU<i)  when 

II  90 

that  termination  occurred  when  max,-  ^(DMaD'1)  =  p(MaUd)  within  four  significant 
figures. 

Termination  conditions  for  supQi  p(MaUnew)  were  adjusted  such  that  the  optimiza- 
tion successfully  halts  when  all  terms  of  the  numerical  gradient  fall  below  10"2.  Tests 
with  varying  termination  conditions  suggest  that  further  reductions  in  the  numerical 
gradient  produce  little  if  any  improvement  in  p(MaUnew). 

Having  discussed  the  three  selection  criteria,  the  flop  count  comparison  of  Table 
8.7  may  be  examined.  The  superiority  of  Method  B  in  terms  of  computing  efficiency 
is  immediately  apparent.  For  the  4  x  4  systems,  an  average  of  172, 166  flops  was 
required  by  Method  A  versus  126,862  flops  for  Method  B:  a  savings  of  26%.  This 
savings  increases  to  42%  for  the  5  x  5  systems  and  63%  for  the  9  x  9  systems.  For 
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design  problems  requiring  exact  values  of  ju,  Method  D  offers  a  savings  of  58%  over 
Method  A  and  39%  over  Method  C  for  9  x  9  systems. 

The  efficiency  of  Method  D  for  finding  the  exact  value  of  y,  seems  to  be  a  direct 
result  of  the  excellent  starting  point  provided  by  first  solving  for  supa.  p(MaUnew). 
In  fact,  for  many  of  the  20  systems  the  flop  count  difference  between  Methods  B 
and  D  simply  reflects  the  calculations  involved  in  computing  ||»  one  time.  The 
resulting  derivative  was  sufficiently  small  to  meet  the  termination  criteria  without 
further  optimization. 

It  is  important  to  emphasize  that  the  criteria  chosen  for  the  flop  count  com- 
parison, especially  the  second  criterion,  actually  skew  the  results  against  the  new 
2(?~  1) -dimensional  search  method.  While  the  flop  count  savings  are  quite  substan- 
tial, the  real  advantages  occur  for  systems  that  do  not  directly  converge  to  ft  from 
the  starting  point.  Convergence  to  some  local  maxima  requires  the  selection  of  a  new 
starting  point  from  which  the  search  for  jjl  may  continue.  For  previous  methods  this 
means  a  search  over  n  —  1  dimensions  so  the  search  area  grows  exponentially  as  sys- 
tem size  increases.  The  new  method,  however,  requires  a  search  over  only  a  2(q  -  1) 
dimensional  space  greatly  increasing  the  likelihood  of  quickly  finding  jx. 

8.9     Algorithm  for  Cusping  Systems 

Certainly  the  numerical  efficiency  of  supaj  p(MaUnew)  offers  advantages  even  if 
the  exact  value  of  \i  must  be  found.  For  the  general  cusping  case  of  q  =  2,  the  ability 
to  generate  and  view  3-D  mesh  plots  as  in  Figure  8.4  offers  a  means  of  partially 
overcoming  the  problem  of  nonconvexity.    A  general  algorithm  for  computing  //  to 
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account  for  cusping  and  noncusping  systems  is  now  presented. 

1.  Compute  mfDa(D  MaD'1).  If  MPDA  is  established  then  ^{DMaD-1)  =  p(Ma) 
and  the  algorithm  is  complete. 

2.  If  MPDA  is  not  established  then  no  stationary  point  was  found  and  the  sys- 
tem is  either  a  "kiss"  or  a  "cusp."  Determine  the  first  q  singular  vectors  of 
(DMaD~l)  where  q  is  the  multiplicity  of  W  and  compute  supa.  p(MaUnew).  If 
MPDA  is  now  established  ("kiss")  or  p{MaUnew)  is  within  some  allowable  dis- 
tance of  a{DMaD-x)  at  the  cusp  then  the  procedure  is  complete  and  upper 
bound  a{DMaD~l)  is  used  for  p{Ma). 

3.  If  p(MaUnew)  falls  below  the  desirable  distance  from  ^{DMaD'1)  then  continue 
with  the  PDA  search  over  ^(DMaD'1)  (Method  D  of  Table  8.7)  to  determine 
a  true  stationary  point. 

4.  If  the  stationary  point  from  Step  3  is  still  not  acceptable  then  either  change 
the  initial  settings  of  the  a,-  and  repeat  Steps  2  and  3  or  generate  a  mesh  of 
p(MaUnew)  over  the  a»  as  discussed  in  Section  8.4.  Then  investigate  all  local 
maxima  from  the  mesh  for  the  global  maximum  of  sup^  p(MaUd). 

While  even  the  mesh  is  not  guaranteed  to  contain  p,  the  advantages  offered  by  the 
optimization  supa_  p(MaUnew)  are  a  tremendous  improvement  over  previous  methods 
of  computing  p  for  cusping  systems.  It  should  also  be  mentioned  that  this  approach 
for  computing  p  also  extends  to  cusping  block  structured  systems.  Equation  5.10 
reveals  the  relationships  between  the  elements  oi  Xu...,Xq  and  Yt,...,Yg  which 
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must  hold  even  when  cusping  occurs.  Similar  relationships  can  be  determined  for  any 
general  block  structured  system  so  forming  Xnew  as  linear  combinations  of  Xu . . . ,  Xq 
proceeds  exactly  as  before.  The  only  difference  between  the  element-by-element  un- 
certainty algorithm  and  that  for  general  block  structured  systems  is  in  the  formulation 
of  Unew  which  must  reflect  the  associated  block  configuration.  Otherwise,  the  same 
reduction  in  optimization  variables  with  corresponding  efficiency  benefits  applies. 


CHAPTER  9 
CONCLUSION 

9.1     Summary 

This  dissertation  presents  new  results  in  the  area  of  stability  robustness  analy- 
sis for  multivariable  systems.  Of  primary  significance  are  the  links  developed  between 
similarity  and  nonsimilarity  scaling  techniques  for  computing  the  structured  singular 
value  (fi).  The  results  of  Chapter  5  show  that  no  more  than  2(tt  -  1)  optimization 
variables  are  required  to  compute  [x  using  either  similarity  or  nonsimilarity  scaling. 
While  the  nonsimilarity  scaling  approach  effectively  addresses  element-by-element 
structured  uncertainties,  the  ability  of  similarity  scaling  to  handle  general  block  struc- 
tured uncertainties  with  no  more  than  2(n  —  1)  free  variables  greatly  simplifies  the 
calculation  of  \i  for  this  important  uncertainty  class. 

For  the  general  block  2x2  problem,  the  elimination  of  redundant  variables  leads 
to  a  guarantee  that  the  solution  of  inf^  a{DbMbD^)  equals  fi(Mb)  regardless  of 
whether  or  not  a  stationary  point  occurs  at  the  "inf."  This  eliminates  the  need  to 
employ  one  of  the  computationally  intensive  lower  bound  optimization  routines  for 
this  class  in  the  event  that  the  maximum  singular  value  repeats. 

In  addition  to  the  reduction  in  optimization  variables,  the  link  between  the  two 
scaling  methods  is  shown  in  Chapter  6  to  include  a  direct  relationship  between  scaling 
matrix  D  of  similarity  scaling  and  scaling  matrices  R  and  L  of  nonsimilarity  scaling. 
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The  value  of  this  direct  relationship  is  shown  to  be  most  evident  when  applied  to 
methods  that  compute  a  lower  bound  for  fj,  as  described  in  Chapter  7.  Since  such 
methods  involve  nonconvex  optimizations,  their  solutions  are  initially  only  candidates 
for  //.  Guaranteeing  that  some  solution  actually  equals  ft  requires  showing  that  the 
candidate  and  the  upper  bound  ofa(bMal)-1)  are  equal.  Computing  the  correspond- 
ing D  for  a  particular  fi  candidate  is  straightforward;  however,  the  actual  singular 
value  decomposition  must  then  be  performed  on  D  and  Ma  matrices  with  dimensions 
as  large  as  n2  x  n2.  By  directly  transforming  D  to  the  corresponding  R  and  L,  the 
upper  bound  may  be  determined  using  the  n  x  n  decomposition  ^{Rr1  ML'1  )a(LPR) 
with  a  substantial  reduction  in  floating  point  operations  for  each  step  in  the  required 
frequency  sweep. 

One  of  the  most  important  advantages  of  computing  the  singular  value  upper 
bound  for  //  is  that  any  local  minimum  must  be  the  global  minimum.  As  long  as 
this  minimum  corresponds  to  a  stationary  point,  the  MPDA  theory  guarantees  that 
this  upper  bound  equals  ft.  However,  for  those  cases  that  do  not  achieve  stationarity 
at  the  minimum  ("cusps"),  ft  is  not  achieved  and,  in  fact,  the  actual  value  of  ft  may 
differ  from  the  upper  bound  by  more  than  15%.  Because  of  the  nonconvexity  involved 
in  lower  bound  calculations  and  the  dangers  of  choosing  an  optimistic  estimate  for 
fi,  the  conservative  upper  bound  is  generally  used  in  place  of  ft  when  cusping  occurs. 
This  is  particularly  true  as  system  size  increases  because  the  lower  bound  solution 
domain  increases  directly  with  n,  greatly  impeding  the  search  process. 

By  applying  the  MPDA  principal,  an  estimate  for  pi  requiring  only  2(q  -  1)  op- 
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timization  variables  (where  q  denotes  the  multiplicity  of  the  cusp)  is  proposed  in 
Chapter  8.  The  largest  deviation  found  to  date  between  this  estimate  and  the  actual 
value  of  ft  is  less  than  0.34%:  well  within  most  real-world  engineering  tolerances. 
Reductions  in  floating  point  operations  of  more  than  60%  for  9  x  9  systems  are  ex- 
perienced versus  previous  methods  and  for  the  most  common  cusping  case  of  q  =  2, 
3 -dimensional  mesh  plots  may  be  used  to  easily  locate  candidates  for  ft.  The  fi 
estimate  holds  for  cusping  systems  with  element-by-element  as  well  as  block  struc- 
tured uncertainties  so  any  system  with  structured  uncertainties  can  benefit  from  this 
method.  Also,  the  estimate  can  serve  as  an  initial  starting  point  in  the  event  that 
ft  must  be  determined  exactly.  Such  a  procedure  employing  the  Principal  Direction 
Alignment  theorem  is  shown  in  Table  8.6  to  find  the  actual  value  of  //  with  only  a 
slight  increase  in  floating  point  operations  over  that  required  for  the  estimate. 

9.2     Future  Directions 

9.2.1  Nonsimilarity  Scaling  with  Block  Structured  Uncertainties 

The  results  of  this  dissertation  serve  to  enhance  many  of  the  important  tools 
used  to  analyze  the  robust  stability  properties  of  multivariable  control  systems.  In 
the  process  of  completing  one  particular  objective,  other  promising  areas  frequently 
appear  as  directions  for  additional  research.  One  of  the  more  intriguing  of  such 
areas  emerging  from  this  work  is  the  possibility  of  adapting  the  nonsimilarity  scaling 
technique  to  work  with  general  block  structured  uncertainties. 

For  element-by-element  structured  uncertainties,  the  advantages  of  nonsimilarity 
scaling  are  quite  clear.    For  example,  using  the  similarity  scaling  approach,  System 
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12  of  Table  5.2  required  922,383  flops  to  compute  fi  using  the  reduced  scaling  struc- 
ture. For  the  same  system,  the  nonsimilarity  scaling  approach  required  only  105,061 
flops:  a  reduction  of  more  than  88%.  Unfortunately,  current  nonsimilarity  scaling 
formulations  do  not  treat  general  block  structured  uncertainties  so  the  tremendous 
advantages  cannot  be  realized  for  this  important  uncertainty  class. 

While  uncertainty  blocks  with  arbitrary  dimension  currently  prevent  the  use  of 
nonsimilarity  scaling  methods,  the  direct  relationship  between  D,  R,  and  L  of  Chapter 
6  may  be  applied  to  at  least  one  specific  block  structure.  This  structure  requires 
only  that  all  the  blocks  be  of  a  common  size  so  that  more  than  the  standard  scalar 
uncertainty  elements  are  allowed.  To  illustrate  the  application  of  nonsimilarity  scaling 
to  this  class  of  uncertainties,  consider  the  system  of  Example  5.4  on  page  61.  A  direct 
calculation  of  R,  and  L  from  D  gives 

R  =  diag[l,    1,   0.833,   0.833],     L  =  diag[l,    1,   0.625,    0.625] 

with  Tj(BrxML-l)a(LPR)  =  a{DbMbD^)  =  16.43  =  p(Mb).  Since  the  direct  trans- 
formation holds  here,  the  problem  could  have  been  solved  by  employing  nonsimilarity 
scaling  from  the  start.  The  ability  of  nonsimilarity  scaling  to  address  even  this  limited 
class  of  block  structured  uncertainties  combined  with  the  correspondence  between  op- 
timization variable  requirements  suggests  that  some  general  block  uncertainty  form 
of  nonsimilarity  scaling  may  eventually  be  possible. 
9.2.2  H°°  and  //-Synthesis  Design  Methods 

Two  recent  advances  in  the  design  of  multivariable  control  systems  are  the 
H°°  and  //-synthesis  methods.  The  H°°  technique  allows  for  the  design  of  controllers 
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that  possess  optimal  performance  in  terms  of  minimizing  the  H°°  norm  of  some  fre- 
quency dependent  transfer  matrix  [18].  On  the  other  hand,  /i-synthesis  extends 
the  optimal  performance  properties  of  the  H°°  method  to  encompass  robust  perfor- 
mance characteristics  as  well  [19].  A  number  of  interesting  design  studies  involving 
the  //-synthesis  approach  have  appeared  including  a  flexible  antenna  structure  and 
an  implementation  of  the  space  shuttle  lateral  axis  flight  control  system  [47,  48].  The 
next  few  paragraphs  give  a  brief  overview  of  these  techniques  along  with  possible 
applications  of  the  results  of  this  dissertation. 

For  transfer  matrix  G,  the  H°°  norm  of  G  may  be  written  as 

\\G{s)\\QO  =  supa(G(s))  (9.1) 

where  s  again  denotes  frequency  dependence    [49].   Figure  9.1  depicts  the  standard 
#°°  system  model  with  feedback  controller  K.  Transfer  matrix  G  can  be  partitioned 


as 

G\i    G\2 
G  = 

G21    G22 

giving  a  set  of  system  equations  of  the  form 


(9.2) 


z  =  Gnw  +  Gi2u,       y  =  G2iw  +  G22U. 

Noting  that  u  =  Ky  allows  this  to  be  rewritten  as 

z  =  [Gn  +  G12K(I  -  G22K)-lG2i]w  (9.3) 

Defining  F  =  [Gn  +  G12K(I  -  G22K)~1G21]  as  the  transfer  matrix  from  external 
input,  w,  to  output  z,  Equation  9.3  can  then  be  written  as 

z  =  Fw. 
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The  design  objective  then  involves  minimizing  the  H°°  norm  of  F.  Depending  the  on 

the  performance  objective(s)  to  be  optimized,  the  choice  of  K  that  best  meets  the 
objective(s)  can  be  chosen  from  the  set  of  all  stabilizing  controllers.  Characteriza- 
tion of  all  proper,  stabilizing  controllers  can  be  obtained  from  the  so-called  "Youla 
parameterization"    [50]  allowing  the  minimization  to  be  formulated  as 


sminKSupa(F(5))  =  smmAJ|F(5)||0O. 


(9.4) 


This  problem  turns  out  to  be  convex  and  a  method  to  directly  solve  for  the  optimal 
controller  has  been  proposed  by  Doyle  et  al    [51]. 


External 
Inputs      w 


Control 

Signals 


4-    Outputs 


Measured 
*     Variables 


Figure  9.1:     H°°  System  Representation. 

The  optimal  performance  offered  by  a  system  satisfying  Equation  9.4  depends 
on  the  adequacy  of  the  plant  model.  Uncertainty  in  this  model  may  lead  to  unpre- 
dictable behavior  that  quickly  loses  optimal  and  even  acceptable  performance  levels. 
However,  by  combining  the  H°°  design  procedure  with  singular  value  scaling  tech- 
niques, robust  performance  and  stability  objectives  may  be  simultaneously  met.  Such 
a  methodology  was  introduced  by  Doyle  under  the  name  "/^-synthesis"    [19].  Rather 
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than  simply  minimizing  a(F(s))  as  a  function  of  frequency  and  stabilizing  controller, 
K,  //—synthesis  allows  for  the  incorporation  of  structured  plant  uncertainties  with 
the  design  objective  of  solving 

min   infllDFD-ML  (9.5) 

stable  K   D    "  "  V        ' 

where  F  represents  the  appropriate  combination  of  F  and  A.  For  fixed  D,  Equation 
9.5  is  a  standard  H°°  problem.  Similarly  for  fixed  K  the  problem  involves  finding 
fi{DFD~l).  Separately,  each  problem  is  convex.  Unfortunately,  the  combined  prob- 
lem loses  this  desirable  convexity  property  so  there  is  no  guarantee  of  ever  achieving 
the  global  minimum  in  terms  of  both  K  and  D. 

While  the  //—synthesis  approach  provides  controllers  with  robust  performance,  it 
requires  a  tremendous  amount  of  computing  power  to  find  even  a  local  minimum  of 
Equation  9.5.  With  the  loss  of  convexity,  computational  requirements  increase  even 
more.  However,  by  applying  the  results  of  this  dissertation,  it  may  be  possible  to 
improve  the  efficiency  of  the  //—synthesis  method. 

For  example,  the  reduction  of  optimization  variables  for  systems  with  both  scalar 
and  block  structured  uncertainties  should  directly  apply  to  that  step  of  the  procedure 
involved  in  infimizing  the  maximum  singular  value  over  D.  Also,  the  direct  relation- 
ship between  similarity  scaling  matrix,  D,  and  nonsimilarity  scaling  matrices  L  and  R 
might  provide  insight  into  a  formulation  of  //—synthesis  involving  only  nonsimilarity 
scaling  with  its  corresponding  reduced  system  size.  Finally,  the  procedure  outlined 
in  Chapter  8  for  estimating  fi  with  a  repeated  maximum  singular  value  should  apply 
to  the  //—synthesis  problem  when  cusping  occurs  so  that  unnecessary  conservatism 
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can  be  reduced  without  the  huge  computing  cost  involved  with  previous  methods  of 
determining  the  lower  bound  for  //. 

Although  the  actual  benefits  of  applying  the  results  of  this  dissertation  to  the 
fi— synthesis  design  approach  are  as  yet  unknown,  any  efficiency  advantage  over  pre- 
vious methods  is  certainly  welcome.  Since  all  but  the  simplest  of  design  problems 
require  tradeoffs  between  performance  requirements,  an  iterative  process  is  normally 
employed  to  achieve  a  satisfactory  performance  compromise.  Efficiency  improvements 
as  low  as  ten  or  twenty  percent  could  therefore  translate  into  several  hours  worth  of 
CPU  savings  as  an  acceptable  performance  trade-off  is  reached. 
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